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train  interactions  is  made  and  features  of  these  types  of  interactions  are  discussed.  In  co 
junction  with  a  review  of  the  existing  literature,  these  results  help  establish  the  type  of 
shock  system  to  be  expected  under  various  operating  conditions  and  the  amount  of  flow  nonuni 
formity  and  unsteadiness  which  can  be  caused. 

The  experimental  investigation  has  been  coordinated  with  a  numerical  study  of  the  shock 
train  phenomenon.  The  computations  use  the  explicit,  time-dependent,  second-order  accurate 
MacCormack  scheme  to  solve  the  mass  averaged  Navier-Stokes  equations.  Turbulence  is  modeled 
via  the  Baldwin-Lomax  algebraic  model  and  the  Wi 1 cox-Rubesi n  two-equation  model.  Test  calcu 
lations  have  been  performed  for  two  flat  plate  equilibrium  turbulent  boundary  layer  flows  an 
one  separated  oblique  shock/turbulent  boundary  layer  interaction.  A  numerical  simulation  of 
the  Mach  1.6  multiple  normal  shock/turbulent  boundary  layer  interaction  is  currently  being 
performed. 
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ABSTRACT 


Multiple  shock  wave/turbulent  boundary  layer  interactions  in  constant  or  nearly 
constant  area  supersonic  duct  flows  occur  in  a  variety  of  devices  including  scramjet  inlets, 
gas  ejectors,  and  supersonic  wind  tunnels.  For  sufficiently  high  duct  exit  pressures,  a 
multiple  shock  wave/turbulent  boundary  layer  interaction  or  shock  train  may  form  in  the 
duct  and  cause  a  highly  nonuniform,  and  possibly  unsteady,  flow  at  the  duct  exit. 

In  this  report,  the  mean  flow  characteristics  of  two  shock  train  interactions,  one 
with  an  initial  Mach  number  of  2.5  the  other  at  Mach  1 .6,  are  investigated  using  spark 
Schlieren  photography,  surface  oil  flow  visualization,  and  mean  wall  pressure 
measurements.  An  LDV  investigation  of  the  Mach  1.6  interaction  is  currently  underway. 
The  experiments  were  performed  in  3  inch  wide  rectangular  ducts  with  nominal  heights  of 
1.5  and  1.4  inches  for  the  Mach  2.5  and  1.6  interactions,  respectively.  The  Mach  2.5 
interaction  was  oblique  and  asymmetric  in  nature.  A  large  separation  occurs  after  the  first 
oblique  shock.  The  top  and  bottom  wall  boundary  layer  separation  has  been  investigated, 
revealing  that  the  shape  of  the  reattachment  lines  and  surface  flow  patterns  for  the  two 
separation  regions  are  quite  different.  This  oblique  shock  flow  pattern  occurs  in  a  neutrally 
stable  fashion  with  each  type  of  opposing  separation  region  alternately  existing  on  either  the 
top  or  bottom  wall  during  the  course  of  a  run.  A  small  scale  unsteadiness  in  the  shock  train 
location,  with  movement  on  the  order  of  a  boundary  layer  thickness,  is  also  observed.  In 
contrast,  the  Mach  1.6  interaction  consists  of  repeated,  symmetric  normal  shocks.  The 
initial,  bifurcated  normal  shock  has  a  small  separation  region  at  its  foot  while  the  following 
weaker  shocks  do  not  separate  the  boundary  layer.  The  number  of  shocks  in  the  shock 
train  and  the  overall  length  of  the  interaction  increase  as  the  boundary  layer  thickens  in  the 
duct.  Only  very  slight  unsteadiness  in  the  shock  train  location  is  observed  at  this  lower 
Mach  number.  A  comparison  of  the  Mach  2.5  and  1.6  shock  train  interactions  is  made  and 
features  of  these  types  of  interactions  are  discussed.  In  conjunction  with  a  review  of  the 
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existing  literature,  these  results  help  establish  the  type  of  shock  system  to  be  expected 
under  various  operating  conditions  and  the  amount  of  flow  nonuniformity  and  unsteadiness 
which  can  be  caused. 

The  experimental  investigation  has  been  coordinated  with  a  numerical  study  of  the 
shock  train  phenomenon.  The  computations  use  the  explicit,  time-dependent,  second-order 
accurate  MacCormack  scheme  to  solve  the  mass  averaged  Navier-Stokes  equations. 
Turbulence  is  modelled  via  the  Baldwin-Lomax  algebraic  model  and  the  Wilcox-Rubesin 
two-equation  model.  Test  calculations  have  been  performed  for  two  flat  plate  equilibrium 
turbulent  boundary  layer  flows  and  one  separated  oblique  shock/turbulent  boundary  layer 
interaction.  A  numerical  simulation  of  the  Mach  1.6  multiple  normal  shock/turbulent 
boundary  layer  interaction  is  currendy  being  performed. 
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I.  INTRODUCTION 


This  report  describes  the  results  of  an  integrated  numerical  and  experimental 
investigation  of  multiple  shock  wave/turbulent  boundary  layer  interactions  in  confined 
rectangular  ducts.  Multiple  shock  wave/turbulent  boundary  layer  interactions,  also  called 
pseudo-shock  or  shock  train  interactions,  are  typically  found  in  confined,  internal,  supersonic 
flows  experiencing  large  adverse  pressure  gradients.  Rather  than  recompressing  through  a 
single  shock,  as  would  be  expected  in  an  inviscid  flow,  the  flow  recompresses  through  an 
extended  shock  system  consisting  of  a  series  of  repeated  normal  or  oblique  shocks. 
Reacceleration  of  the  flow  following  each  shock  is  caused  by  adjustments  in  the  wall  boundary 
layers.  The  objective  of  this  study  is  to  investigate  the  shock  train  phenomenon,  both 
experimentally  and  numerically,  with  the  goal  of  understanding  the  detailed  turbulent  and  mean 
flow  mechanisms  occurring  in  the  shock  train.  Additionally,  the  ability  of  existing  numerical 
techniques  to  compute  such  flows  is  evaluated. 

Multiple  shock  wave/turbulent  boundary  layer  interactions  occur  in  a  variety  of  devices 
of  technological  importance.  Examples  include  supersonic  gas  ejectors,  supersonic  wind 
tunnel  diffusers,  and  scramjet  (supersonic  combustion  ramjet)  inlets.  Ejectors  are  used  to  both 
pump  and  mix  fluids  and  are  especially  suited  to  applications  requiring  low  maintenance  or 
involving  corrosive  fluids.  A  supersonic  ejector  uses  a  high  stagnation  pressure  supersonic 
stream  to  entrain  a  lower  stagnation  pressure  secondary  stream.  The  two  streams  mix  and 
recompress  in  a  mixing  duct.  Under  proper  conditions,  both  streams  attain  supersonic 
velocities  in  the  duct,  then  decelerate  through  a  shock  train  system  to  subsonic  exit  velocities. 
Both  the  mixing  effect  of  the  shock  train  as  well  as  its  pressure  recovery  are  important  in  this 
application.  In  supersonic  wind  tunnel  diffusers,  shock  train  systems  are  often  found  and  are 
usually  the  dominant  source  of  losses  in  the  diffuser.  Consequently,  the  power  required  to 
operate  the  wind  tunnel  is  directly  related  to  the  losses  occurring  in  the  shock  train. 
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The  main  impetus  for  the  current  research  is  the  development  of  scramjet  propulsion 
systems.  Recently,  a  heightened  interest  in  scramjet  propulsion  has  been  shown.  The  U.S. 
Navy  is  developing  defensive  missiles  capable  of  maintaining  high  supersonic  speeds  for 
ranges  exceeding  several  hundred  miles.  Ramjet  or  scramjet  air-breathing  propulsion  systems 
will  probably  be  required  to  achieve  the  necessary  range  under  the  size  and  weight  constraints 
imposed  by  current  launch  platforms  and  launchers.  Much  attention  has  also  been  focused  on 
the  development  of  a  new  hypersonic  aircraft,  called  the  National  Aero-Space  Plane  or  X-30, 
w'hich  will  employ  scramjet  propulsion1.  The  development  of  such  an  aircraft  is  expected  to 
depend  heavily  on  computational  fluid  dynamic  (CFD)  methods  to  supplement  conventional 
wind  tunnel  testing,  assuming  the  reliability  of  the  computational  approach  can  be 
demonstrated2.  This  increased  dependence  on  numerical  predictions  is  due  to  the  scarcity  of 
hypersonic  test  facilities,  difficulties  in  performing  hypersonic  measurements,  and  the  lower 
relative  cost  of  numerical  predictions.  The  research  reported  herein  will  help  develop  a  bener 
understanding  of  the  inlet  flow,  including  the  detailed  turbulent  transport  and  flow 
reacceleration  phenomena  which  can  occur  and  will  help  qualify  computational  techniques  for 
predicting  such  flows. 

Traditionally,  less  than  optimal  scramjet  and  ramjet  inlet  designs  have  been  used  due 
to  the  limited  knowledge  of  the  inlet  flowfield.  Ramjet  inlet  flowfields  are  characterized  by 
complicated  shock  structures,  rapidly  growing  boundary  layers,  large  separated  regions,  and 
undesirable,  self-excited,  low  frequency  oscillations.  Ramjet  inlet  flows  have  also  shown 
strong  coupling  with  unsteady  combustor  pressures  resulting  in  large  amplitude  oscillations. 
These  problems  have  led  to  a  series  of  investigations,  both  numerical  and  experimental,  into  the 
behavior  of  ramjet  flowfields.  Some  of  the  more  resent  results  are  those  of  Hsieh,  et  alA 
Hsieh,  et  al.4,  Bogar5,  Bogar,  et  al.6,  Sajben,  et  al.7,  and  Talcott  and  Kumar8.  Scramjet  inlet 
Hows  have  received  considerably  less  attention  relative  to  those  of  ramjet  inlets.  While  the  two 
types  of  inlets  perform  similar  functions,  the  ramjet  inlet's  geometry  and  internal  flow  are 
significantly  different  than  the  scramjet's.  Scramjet  inlet  flowfields  are  also  characterized  by 
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complicated  shock  patterns,  large  separation  regions,  and  thick  boundary  layers  compared  to 
the  duct  dimensions.  For  sufficiently  high  combustor  pressures,  a  shock  train  may  form  in  the 
diffuser  portion  of  the  inlet,  thereby  adversely  affecting  the  pressure  recovery  and  causing  a 
highly  nonuniform  flow  at  the  diffuser  exit.  The  unsteady  nature  of  the  scramjet  inlet  flow,  the 
inlet  normal  and  oblique  shock  train  phenomena,  the  susceptibility  of  the  scramjet  inlet  to  self- 
excited  and  combustor-induced  oscillations,  and  the  unsteady  inlet  buzz  phenomenon  are 
aspects  of  the  scramjet  inlet  flow  which  require  further  study.  In  particular,  Waltrup9  identified 
the  shock  wave/boundary  layer  interactions  occurring  in  scramjet  inlets  as  an  area  requiring 
further  research.  He  states  that  detailed  measurements  of  the  turbulent  transport  and  dissipation 
mechanisms  in  the  near  wall  region  for  a  variety  of  initial  conditions  are  needed.  According  to 
Waltrup,  "it  is  these  that  are  of  interest  in  hypersonic  inlets,  especially  in  internally  ducted 
supersonic  flows  with  compression  and  expansion  waves  so  that  the  viscous  total  pressure 
losses  as  well  as  regions  of  separation  can  be  modeled  and  predicted  with  confidence." 

As  stated  above,  the  objective  of  this  research  is  to  contribute  to  a  better  understanding 
of  the  complicated  flow  mechanisms  involved  in  the  shock  train  phenomenon  and  to  evaluate 
numerical  predictive  techniques  for  such  flows.  A  combined  numerical  and  experimental 
investigation  has  been  performed.  The  numerical  portion  of  the  work  was  undertaken  both  to 
evaluate  the  ability  of  existing  numerical  techniques  for  calculating  such  flows  and  to  assist  in 
the  planning,  execution,  and  evaluation  of  the  experimental  work.  These  calculations  employ 
the  widely  used  explicit  MacCormack  scheme  to  integrate  the  time-dependent,  mass-averaged 
Navier-Stokes  equations  along  with  a  two-equation  model  of  turbulence  and  an  algebraic 
turbulence  model.  This  approach  is  typical  of  those  currently  being  used  for  scramjet  engine 
analysis  as  discused  by  White,  et  al.1®,  although  most  previous  work  has  used  only  the  simpler 
algebraic  turbulence  models  and  has  focused  on  the  gross  flowfield  characteristics.  The  current 
experimental  investigation  has  been  performed  to  improve  the  understanding  of  the  flow 
mechanisms  occurring  in  the  shock  train  interaction  and  to  add  to  the  relatively  sparse  data  base 
on  such  interactions.  Such  a  data  base  is  needed  both  to  assist  in  the  understanding  and 
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analysis  of  shock  train  interactions  and  to  evaluate  and  improve  numerical  techniques  for 
calculating  such  flows.  The  experimental  measurements  were  taken  in  a  small  scale  planar 
two-dimensional  wind  tunnel  with  incoming  Mach  numbers  of  2.45  and  1.6.  Experimental 
measurement  techniques  included  spark  Schlieren  photography,  oil  streak  surface  flow 
visualization,  mean  wall  static  pressure  measurements,  and  two-component,  coincident  laser 
Doppler  velocimeter  (LDV)  measurements  of  the  mean  and  turbulent  velocity  fields.  The 
experimental  techniques  employed  are  similar  to  those  being  used  by  Yanta11  to  investigate  a 
planar  tw'o-dimensional  scramjet  inlet  flow.  This  inlet  consists  of  a  10  degree  precompression 
ramp,  followed  by  an  inward  turning  scoop  feeding  a  nearly  constant  area  supersonic  diffuser 
of  rectangular  cross  section.  Waltrup9  listed  Yanta’s  work  as  the  only  other  relevant 
investigation  of  a  diffuser  flow  of  this  type  in  which  experimental  measurements,  of  sufficient 
detail  to  study  the  turbulent  transport  as  well  as  mean  flow  phenomena  occurring,  are  being 
taken.  In  the  current  work,  the  flow  non-uniformities  caused  by  the  inward  turning  scoop 
portion  of  a  scramjet  type  inlet  have  been  avoided  and  the  flow  mechanisms  in  the  shock  train 
system  are  investigated  for  a  uniform  incoming  supersonic  ducted  flow,  with  an  equilibrium 
turbulent  boundary  layer. 

Several  studies12’13, 14  have  shown  that  three-dimensional  effects  are  present  in 
nominally  two-dimensional  planar  oblique  shock/boundary  layer  interactions  due  to  the 
interaction  of  the  oblique  shock  with  the  side  wall  boundary  layers.  Chriss,  et  al.15  have  also 
shown  that  some  three-dimensionality  is  present  in  a  multiple  normal  shock  interaction  in  a 
square  duct.  While  axisymmetric  geometries  provide  a  more  nearly  two-dimensional  flow, 
measurements  in  these  geometries  are  difficult.  Optical  techniques  such  as  Schlieren  or  LDV 
are  complicated  by  the  curvature  of  the  duct  walls.  Physical  probes,  such  as  pressure  or  hot 
wire  probes,  can  be  used,  but  their  presence  causes  undesirable  flow  disturbances.  In  this 
study,  a  planar  two-dimensional  geometry  was  selected  to  facilitate  LDV  and  Schlieren  studies 
of  the  interaction.  Furthermore,  applications  such  as  the  National  Aerospace  Plane  will  involve 
planar  geometries. 
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In  the  remainder  of  this  report,  the  shock  train  phenomenon  will  be  described  in  more 
detail.  The  literature  review  describes  both  confined  and  unconfined  shock  wave/boundary 
layer  interactions  illustrating  the  major  differences  between  these  two  types  of  flows.  The 
major  features  of  the  confined,  shock  train  interaction  are  introduced  and  previous  experimental 
and  numerical  investigations  of  this  phenomenon  are  discussed.  Details  of  the  experimental 
and  numerical  approaches  are  then  given.  Finally,  a  discussion  of  the  experimental  and 
numerical  results  is  presented  and  some  concluding  remarks  and  recommendations  are  made. 


II.  LITERATURE  REVIEW 


In  this  section  the  relevant  literature  pertaining  to  multiple  shock  wave/turbulent 
boundary  layer  or  shock  train  interactions  in  constant  area  or  nearly  constant  area  rectangular 
ducts  is  reviewed.  An  introduction  to  the  nomenclature  and  basic  physical  concepts  relating  to 
this  type  of  interaction  is  incorporated  into  the  literature  review.  Before  proceeding  further,  a 
distinction  between  confined  and  unconfined  shock  wave/turbulent  boundary  layer  interactions 
is  made.  The  confined  nature  of  an  internal  ducted  supersonic  flow  is  necessary  for  the  shocr. 
train  phenomenon  to  occur.  With  external  shock  wave/boundary  layer  interactions,  the  flow 
downstream  of  the  interaction  is  unconfined,  and  the  multiple  shock  phenomenon, 
characterized  by  flow  reacceleration  following  each  shock,  is  not  normally  observed.  Thus, 
the  shock  wave/boundary  layer  interaction  for  confined  flows  is  fundamentally  different  than 
for  external,  unconfined  flows.  In  this  work  a  confined  interaction  is  defined  as  an  internal 
interaction  in  which  the  confinement  effect  of  the  duct  walls  noticeably  affects  the  interaction. 
Unconfined  interactions  are  defined  as  any  interaction,  external  or  internal,  in  which 
confinement  effects  are  negligible. 

A  large  number  of  previous  studies  have  focused  on  the  shock  wave/turbulent 
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boundary  layer  interaction.  Fortunately,  several  review  articles  are  available  which  have 
summarized  the  results  of  these  papers.  Two  review  articles  are  of  particular  interest:  Green16 
and  Adamson  and  Messiter17.  While  no  effort  will  be  made  to  duplicate  the  contents  of  these 

! 

two  reviews,  several  papers  deserve  individual  discussion  due  to  their  relevance  to  the  work 
reported  herein.  Additionally,  important  contributions  to  this  area  of  research  have  been  made 
since  Adamson  and  Messiter's17  article. 

t 

A.  EXPERIMENTAL  WORK,  UNCONFINED  INTERACTION 

In  this  section,  only  papers  treating  the  unconfined,  nominally  steady,  normal 
shock/turbulent  boundary  layer  interaction  will  be  reviewed,  as  this  unconfmed  interaction 


bears  some  resemblance  to  the  confined  shock  train  interaction  under  consideration.  A 
summary-  of  the  relevant  unconfined  studies  is  given  in  Table  1.  Several  researchers  have 
considered  this  problem  including  Seddon18,  Vidal,  et  al.19,  Kooi211,  East21,  and  Delery22 
with  both  East  and  Delery  making  laser  Doppler  velocimeter  (LDV)  measurements  of  the  flow . 
These  studies  were  all  performed  in  wind  tunnels  for  which  confinement  effects  were  small. 
The  flows  studied  by  East  and  Delery  exhibited  small  flow  confinement  effects  and  one  of 
Seddon's  configurations  also  showed  similar  effects.  However,  these  three  studies  18, 2 1.22 
appear  to  be  best  described  as  unconfined  interaction  studies. 

East21  and  Delery22  both  observed  unseparated  normal  shock/turbulent  boundary  layer 
interactions  with  a  Mach  number  at  the  start  of  the  interaction,  Mu,  of  1.3.  At  this  low  Mach 
number,  the  pressure  rise  across  the  shock  is  too  small  to  cause  the  boundary  layer  to  separate. 
As  shown  in  Figure  1,  the  subsonic  portion  of  the  unseparated  boundary  layer  thickens 
upstream  of  the  shock  causing  compression  waves  to  emanate  from  the  sonic  line  which 
eventually  coalesce  forming  a  foot  for  the  outer  shock.  The  outer  shock  is  slightly  curved, 
being  a  strong  oblique  shock  (not  strictly  normal)  with  the  flow  following  the  interaction  being 
totally  subsonic.  The  wall  static  pressure  increases  in  a  continuous  fashion,  while  the  static 
pressure  in  the  core  flow  increases  discontinuously,  reflecting  the  mixed  subsonic/supersonic 
nature  of  the  interaction. 

At  higher  Mach  numbers,  the  inertial  forces  in  the  subsonic  layer  near  the  wall  are  not 
strong  enough  to  negotiate  the  increased  pressure  rise  and  the  boundary  layer  will  separate  at 
the  foot  of  the  shock,  thereby  increasing  the  streamwise  extent  of  the  pressure  rise.  Seddon18, 
Vidal  et  al.19,  Kooi2^,  East21,  and  Delery22  all  found  separated  normal  shock/turbulent 
boundary  layer  interactions  at  Mach  numbers  greater  than  Mu  =  1.4  and  various  Reynolds 
numbers.  The  presence  of  a  free  interaction  at  the  separation  point  has  been  observed 
experimentally,  i.e.  the  flow  upstream  of  and  slightly  downstream  of  the  separation  does  not 
depend  on  the  details  of  the  downstream  pressure  rise. 


A  schematic  of  the  separated,  unconstrained  interaction  is  shown  in  Figure  2.  The 
increasing  thickness  of  the  subsonic  layer  in  the  region  near  the  separation  point  causes  weak 
oblique  compression  waves  to  propagate  into  the  outer  flow.  These  waves  coalesce  to  form  a 
weak  oblique  shock  (leading  shock)  which  eventually  intersects  the  strong,  nearly  normal 
shock  in  the  outer  supersonic  flow.  At  this  intersection,  termed  the  bifurcation  point,  a  second 
weak  oblique  shock  (trailing  shock)  is  generated  which  propagates  back  toward  the  w'all.  A 
slip  line  is  also  generated  at  the  bifurcation  point  and  extends  downstream,  static  pressure  and 
flow  direction  being  matched  across  the  slip  line.  The  second  oblique  shock  is  required  to 
satisfy  compatibility  of  flow  direction,  while  the  slip  line  is  indicative  of  a  mismatch  in  velocity 
magnitude  caused  by  the  differing  losses  through  the  outer  and  inner  shocks.  A  region  of 
supersonic  flow  which  isentropically  decelerates  to  subsonic  speed  may  exist  in  the  otherwise 
subsonic  flow  below  the  slip  line.  In  experiments  involving  a  Mach  1.47  separated  normal 
shock  interaction,  Seddon18  observed  such  a  region,  which  he  called  a  "supersonic  tongue." 
However,  in  a  similar  separated  normal  shock  interaction  at  Mach  1.4,  Kooi20  did  not  observe 
this  phenomenon,  indicating  a  Mach  number  dependence.  This  dependence  was  confirmed  in  a 
series  of  three  normal  shock/boundary  layer  interaction  experiments  where  East21  found  all 
subsonic  flow  behind  a  Mach  1.3  interaction,  a  region  of  sonic  flow  behind  a  Mach  1.4 
interaction,  and  a  supersonic  tongue  at  the  boundary  layer  edge  behind  a  Mach  1.54  interaction. 
A  second  Mach  number  dependence  is  related  to  the  strength  of  the  nearly  normal  shock  As 
the  Mach  number  increases,  the  pressure  rise  increases,  causing  the  separation  region  to  grow 
in  size,  effectively  pushing  the  bifurcation  point  further  from  the  wall.  The  freestream 
Reynolds  number  also  has  an  effect  on  the  interaction  as  shown,  for  example,  by  Vidal,  et 
al.19.  For  a  fixed  pressure  rise,  the  length  of  the  separation  region  and  height  of  the  bifurcation 
point  above  the  wall  increases  as  Reynolds  number  based  on  undisturbed  boundary  layer 
thickness.  Regu,  decreases. 

Some  general  observations  about  these  unconfined,  normal  shock  wave/turbulent 
boundary  layer  studies  follow.  Throughout  the  interaction  region,  large  variations  in  the  static 


and  stagnation  pressure  are  observed  across  the  boundary  layer,  normal  to  the  wall,  and 
following  the  interaction  the  boundary  layer  takes  an  appreciable  distance  to  fully  recover  to 
equilibrium  conditions.  For  example,  Seddon1^  found  the  boundary  layer  following  a 
separated  normal  shock  interaction  to  be  still  recovering  after  up  to  50  undisturbed  boundary 
layer  thicknesses  downstream  of  the  start  of  the  interaction.  Vidal,  et  al.19,  using  high  speed 
Schlieren  movies,  detected  an  unsteady  shock  structure  at  the  bifurcation  point.  Unsteady 
shock  structures  were  also  observed  by  East21  and  by  Delery22.  Both  investigators  found  that 
the  use  of  an  adjustable  second  throat  to  choke  the  subsonic  flow  following  the  interaction 
substantially  improved  the  steadiness  of  the  shock,  presumably  because  such  a  device  isolates 
the  interaction  from  downstream  pressure  fluctuations.  East21  also  used  conditional  sampling 
in  his  LDV  measurements  to  avoid  measuring  "false  turbulence"  caused  by  shock  motion. 

B.  EXPERIMENTAL  WORK,  CONFINED  INTERACTION 

Despite  many  similarities,  major  differences  exist  between  unconfined  and  confined 
shock  wave/boundary  layer  interactions.  Flow  confinement,  as  characterized  by  the  ratio  of  the 
undisturbed  boundary  layer  thickness  to  the  radius  of  the  duct  (for  axisymmetric  geometries), 
§u/r,  can  have  significant  effects  on  the  interaction,  leading  to  three  types  of  shock  systems23:  a 
single  normal  shock  for  small  8u/r,  a  series  of  nearly  normal  shocks  for  moderate  5u/r,  and  a 
senes  of  strong  oblique  shocks  for  large  Sj/r.  In  a  planar  geometry,  the  confinement  parameter 
is  defined  as  Sy/h,  where  h  is  the  duct  half  height  in  the  direction  perpendicular  to  the  surface  at 
which  the  interaction  of  interest  occurs  (normally  taken  as  the  shorter  dimension  in  a 
rectangular  duct).  The  Mach  number  in  the  core  flow  at  the  start  of  the  inteiaction,  Mu.  and  the 
confinement  parameter  are  the  primary  parameters  affecting  the  confined  interaction  as  shown 
by  \lerkli2J  and  Mateer  and  Viegas2^.  Reynolds  number  was  shown  to  have  a  less 
pronounced  influence  on  the  confined  interaction  as  compared  to  the  unconfined  interaction. 
Experimental  studies  of  the  confined  shock  wave/boundary  layer  interaction  have  been 
performed  in  both  planar  two-dimensional  and  axisymmetric  geometries,  with  varying  degrees 


of  detail.  In  reviewing  the  experimental  work,  the  various  types  of  measurements  will  be 
discussed  separately,  beginning  with  the  Schlieren  flow  visualization,  overall  pressure 
recovers’,  and  surface  pressure  rise  studies,  and  followed  by  the  mean  flow  and  turbulence 
studies.  A  summary  of  the  confined  experimental  studies  is  given  in  Table  2. 

Lustwerk23,  McLafferty,  et  al.2^,  and  Fejer,  et  al.27  have  all  shown  still  Schlieren 
photographs  of  the  various  shock  structures.  While  these  photographs  clearly  show  the 
instantaneous  structure  of  the  three  types  of  shock  systems  (single  normal  shock,  multiple 
normal  shocks,  or  multiple  oblique  shocks),  high  speed  Schlieren  movies  of  the  interactions 
were  not  used  to  investigate  the  steadiness  of  the  interaction.  McLafferty,  et  al.26  also 
presented  high  speed  measurements  of  fluctuating  wall  pressures  in  a  constant  area 
axisymmetric  duct.  Pressure  fluctuations  with  frequencies  of  240  to  300  Hz  were  found  in  the 
neighborhood  of  the  shock  interaction.  These  fluctuations  were  present  at  Mach  numbers  of 
2.51  and  2.20,  indicating  the  presence  of  some  shock  unsteadiness  in  constant  area,  circular 
ducts.  Resonance  frequencies  of  about  300  Hz  were  also  calculated  for  the  subsonic  duct  flow 
following  the  shock.  Lustwerk23  noted  that  diverging  the  duct  by  50  minutes  (full  angle) 
helped  to  stabilize  the  shock  system  in  a  rectangular  duct. 

A  series  of  papers  from  by  Ikui,  et  al.28.29.30  gives  an  enlightening  overview  of  the 
shock  train  phenomenon.  Schlieren  photographs  along  with  wall  pressure  measurements  were 
used  to  investigate  the  shock  train  for  Mach  numbers  ranging  from  1.33  to  2. 79. 28  At  Mu  = 
1.33  a  single  normal  shock  was  observed.  Mach  numbers  from  1.37  to  1.60  yielded  multiple 
normal  shocks,  while  Mach  numbers  from  1.86  to  2.79  caused  multiple  oblique  shocks. 
When  the  Mach  number  exceeded  1.60  considerable  flow  asymmetry  was  observed.  Shock 
wave  motion  was  observed  at  all  these  Mach  numbers  and  was  investigated  using  high  speed 
Schlieren  movies  and  high  speed  wall  pressure  measurements.29  The  amplitude  of  the  shock 
oscillation  increased  with  increasing  Mach  number.  They  theorized  that  the  shock  motion  is 
triggered  by  high  frequency  pressure  fluctuations  in  the  incoming  boundary  layer  while  the 
frequencies  of  the  shock  motion  are  related  to  the  resonance  frequencies  of  the  downstream 
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subsonic  duct  flow.  Reference  30  by  Ikui,  et  al.  considered  ducts  with  relatively  large 
divergence  and  is  not,  therefore,  directly  applicable  to  the  current  study. 

The  early  investigators  were  mainly  concerned  with  the  overall  pressure  recovery 
characteristics  of  constant  area  supersonic  diffusers.  Johnson  and  Wu31  compiled  a  summary 
of  these  early  diffuser  pressure  recovery  results.  Typically,  the  static  wall  pressure  was  found 
to  rise  in  a  continuous  fashion  over  a  distance  of  8  to  12  duct  diameters  with  the  pressure 
recovery  slightly  lower  than  that  occurring  for  a  single  normal  shock  at  the  same  undisturbed 
Mach  number.  In  more  recent  results,  Waltrup  and  Billig32,  Merkli24,  Mateer  et  al.33,  Mateer 
and  Viegas25,  Om,  et  al.34,  and  Om  and  Childs35  have  studied  various  parametric  trends  in  the 
pressure  rise.  Some  difficulty  arises  when  comparing  these  investigations  since  the  data  is 
presented  in  various  formats  and  the  parameters  are  not  always  non-dimensionalized  in  the 
same  way.  The  start  of  the  interaction  is  typically  defined  as  the  point  at  which  the  wall 
pressure  begins  to  rise,  and  flow  parameters  at  this  location  are  referred  to  as  undisturbed 
values.  In  the  more  recent  papers,  the  wall  pressure  is  normalized  by  the  undisturbed  pressure. 
Pu,  or  the  upstream  stagnation  pressure,  Pq,  and  the  axial  distance  is  measured  from  the  start  of 
the  interaction  and  normalized  by  the  undisturbed  boundary  layer  thickness,  5U.  Several  trends 
are  noticed  in  the  wall  pressure  distribution  when  plotted  in  such  a  manner.  An  initial  steep 
pressure  rise  is  followed  by  a  region  of  nearly  constant  slope  in  the  pressure  distribution.  The 
shape  of  the  wall  pressure  profile  displays  a  strong  dependence  on  Mach  number,  My,  and 
confinement  parameter,  6u/r,  and  a  weaker  dependence  on  Reynolds  number.  The  slope  in  the 
initial  pressure  rise  increases  slightly  with  increasing  unit  Reynolds  number.  Re/m.  This  result 
seems  to  be  valid  for  both  single  and  multiple  shock  systems  (possibly  because  of  a  free 
interaction  effect  at  the  start  of  the  interaction).  Conflicting  Mach  number  dependencies  were 
found.  For  the  single  shock,  Om,  et  al.34  found  the  slope  of  the  initial  pressure  rise  tends  to 
increase  with  increasing  Mach  number  and  constant  unit  Reynolds  number.  Re/m.  while 
Mateer  and  Viegas25  found  the  initial  slope  to  decrease  with  increasing  Mach  number  and 
constant  Reynolds  number  based  on  downstream  distance,  Res.  The  Reynolds  number  based 


on  undisturbed  boundary  layer  thickness,  Resu,  was  not  held  constant  in  either  experiment. 
Om,  et  al.34  noted  that  Re§u  decreased  at  a  faster  rate,  as  Mach  number  decreased,  in  Mateer 
and  Viegas'  25  series  of  experiments  than  in  their  own,  and  suggest  that  this  is  the  cause  of  the 
conflicting  Mach  number  dependencies.  This  observation  implies  that  the  proper  Reynolds 
number  to  use  in  characterizing  the  flow  is  one  based  on  an  undisturbed  boundary  layer 
thickness.  One  consistent  Mach  number  effect  is  that  the  length  of  the  interaction  increases 
with  increasing  Mach  number.  Considering  the  effect  of  confinement,  one  observes  that  the 
pressure  recovery  decreases  with  increasing  §u/r  and  fixed  Mach  number  and  unit  Reynolds 
number35.  This  effect  may  be  due  to  increased  losses  as  the  shock  system  becomes  a  multi¬ 
shock  structure  and  to  a  greater  increase  in  the  displacement  thickness  for  the  multiple  shocks 
as  opposed  to  the  single  shock  system.  The  increase  in  displacement  thickness  effectively 
accelerates  the  subsonic  flow  after  the  shock  system  causing  a  reduction  in  static  pressure. 

A  description  of  the  local  mean  and  fluctuating  flow  characteristics  will  now  be  given, 
with  the  single  shock  interaction  and  multiple  shock  interaction  being  considered  separately. 
The  confined,  separated  single  normal  shock  interaction  shown  in  Figure  3  is  fairly  similar  to 
the  unconfined,  separated  normal  shock  interaction  (Figure  2).  Mateer  and  Viegas25  showed 
that  the  tendency  toward  separation  increased  as  5u/r  decreased,  i.e.  flow  confinement  delays 
separation.  As  was  found  for  the  unconfined  case,  the  height  of  the  bifurcation  above  the  wall 
and  the  tendency  toward  separation  increases  as  Mach  number,  My,  increases  and  as  unit 
Reynolds  number.  Re/m,  decreases.  Mateer  et  al.33  showed  that  the  length  of  the  separation 
scales  directly  with  the  undisturbed  boundary  layer  thickness.  Their  measurements,  made  with 
the  use  of  a  fast  responding,  embedded  hot  wire  separation  detector,  indicated  some 
unsteadiness  in  the  separation  region.  In  this  study,  surveys  of  the  mean  and  fluctuating  flow 
were  also  made  using  pitot  pressure  and  hot  wire  probes.  Upstream  of  the  shock,  a  peak  in  the 
axial  mass  flux  fluctuations  was  observed  in  the  boundary  layer,  away  from  the  wall.  An  order 
of  magnitude  increase  in  the  mass  flux  fluctuations  was  observed  going  through  the  interaction 
followed  by  a  relaxation  and  diffusion  of  the  fluctuating  quantities  downstream  of  the 
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interaction  leading  to  lower,  nearly  constant  levels  through  the  boundary  layer  out  into  the 
mean  flow.  At  transonic  speeds,  the  interpretation  of  hot  wire  data  is  difficult  due  to  the 
presence  of  shock  waves  on  the  probe.  As  a  result,  Matecr  et  al.33  commented  that  the 
accuracy  of  the  hot  wire  data  may  be  poor  for  the  range  of  local  Mach  numbers,  1 .0  <  M  < 
1.25.  They  also  noted  the  presence  of  a  supersonic  tongue  at  an  initial  Mach  number  of  1.5  but 
did  not  mention  the  presence  of  a  slip  line  in  the  flow.  Om,  et  al.34  made  a  detailed  study  of 
the  single  shock  interaction  in  an  axisymmetric  duct  using  pitot  and  static  pressure  probes  to 
map  out  the  mean  flow,  and  an  alcohol  injection  technique  to  detect  separation.  They  found  no 
separation  at  a  Mach  number  of  1.28,  incipient  separation  at  Mu  =  1.37,  and  separated  flow  at 
Mu  =  1.48.  Their  measurements  also  showed  an  embedded  supersonic  tongue  for  Mu  =  1.28, 
1.37,  and  1.48  which  increased  with  size  as  Mach  number  increased.  The  shape  of  their 
supersonic  tongue  was  different  than  that  found  in  external  interactions 1 8,20,2 1  jn  that  the 
thickness  of  the  tongue  increased  with  downstream  distance.  Om,  et  al.34  observed  that  the 
confinement  effect  tended  to  increase  the  size  of  the  supersonic  tongue.  By  increasing  the  unit 
Reynolds  number,  the  confinement  effect,  and  correspondingly  the  size  of  the  supersonic- 
tongue,  were  reduced.  A  slip  line  was  detected  at  Mu  =  1.48,  but  was  too  weak  to  be 
measured  at  Mach  numbers  of  1.37  and  1.28.  In  one  of  the  few  detailed  studies  performed  at 
higher  Mach  numbers,  Cuffel  and  Back36  investigated  a  Mach  2.92  interaction  using  pressure 
and  temperature  probes.  In  spite  of  relatively  high  confinement  levels,  only  a  single,  separated 
normal  shock  was  reported.  Wall  cooling  was  present  in  this  experiment  and  may  have  had 
sufficient  effect  to  change  the  trends  observed  in  other  adiabatic  wall  experiments. 

Using  the  same  test  section,  Om,  et  al.  34  found  single  normal  shocks  with  Mu  =  1.48. 
Re/m  =  4.92x  1()6,  and  5u/r  =  0.081  while  Om  and  Childs  33  found  a  repeated  nearly  normal 
shock  pattern  at  essentially  the  same  Mach  number  and  Reynolds  number.  Mu  =  1.49.  Re/m  = 
4.90.x  106.  and  higher  confinement  levels.  SJr  =  0.198.  Om  and  Childs35  investigated  this 
multiple  shock  interaction  making  detailed  pitot  and  static  pressure  surveys.  A  qualitative 
sketch  of  a  two  shock  system  is  shown  in  Figure  4.  The  repeated,  nearly  normal  shock 
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interaction  investigated  by  Om  and  Childs35  consisted  of  several  shocks  whereas  the  sketch  in 
Figure  4  shows  only  two  shocks.  In  their  study,  the  boundary  layer  displacement  thickness 
was  found  to  increase  rapidly  through  the  First  shock,  then  decrease  slightly  before  increasing 
again  across  the  second  shock.  This  pattern  was  repeated  through  the  whole  interaction.  The 
momentum  thickness  distribution  exhibited  similar  behavior.  A  slip  line  was  generated  at  the 
bifurcation  point  of  the  first  shock  and  persists  through  the  following,  weaker  shocks  which  do 
not  appear  to  generate  slip  lines.  Following  the  interaction,  the  pitot  pressure  profiles  remain 
nonuniform  even  after  the  static  pressure  profiles  have  become  uniform.  The  flow  just  outside 
the  boundary  layer  stayed  supersonic  through  the  whole  interaction,  while  the  flow  in  the  core 
region  went  subsonic  after  each  shock,  then  reaccelerated  to  supersonic  speed  before  the  next 
shock,  the  distance  between  successive  shocks  decreasing  with  downstream  distance  from  the 
start  of  the  interaction.  At  the  end  of  the  interaction,  the  flow  is  still  mixed  supersonic/subsonic 
in  nature,  with  the  supersonic  portion  of  the  flow  following  the  last  shock  isentropically 
decelerating  to  subsonic  velocities.  Recently,  an  investigation  of  the  three-dimensional  nature 
of  the  shock  wave/turbulent  boundary  layer  interaction  in  a  square  duct  was  made  by  Chriss,  et 
al.15  using  a  one-component  LDV  system  to  make  non-coincident  measurements  of  the  three 
velocity  components.  At  Mach  1 .3  an  unseparated  single  shock  system  was  found  that  was 
largely  two-dimensional  with  the  notable  exception  of  a  supersonic  tongue  near  the  comer.  A 
separated,  multiple  normal  shock  system  was  found  at  Mach  1.6  that  included  a  higher  level  of 
three-dimensionality  especially  near  the  comers  where  the  size  of  the  supersonic  tongue  was 
greater  than  at  mid-span.  The  extent  of  flow  separation  was  also  greatest  in  the  comers. 
Considerable  high  frequency  shock  motion  was  observed  and  caused  distinctly  bimodal 
velocity  histograms  in  the  LDV  measurements  near  the  shock.  Schlieren  photographs  revealed 
what  appears  to  be  a  weak  reflected  shock  generated  where  the  trailing  leg  of  the  initial, 
bifurcated  shock  intersects  the  boundary  layer.  Measurements  were  not  made  in  the  near  wall 
portion  of  the  boundary  layer,  and  boundary  layer  thickness  parameters  were  not  reported  nor 
was  any  turbulence  informatton  apparently  obtained.  Data  for  the  spanwise  component  ot 
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velocity  were  not  included  in  this  paper  either  but  are  expected  to  be  released  at  a  later  date.  In 
the  only  other  detailed  study  of  the  mean  flow  in  a  confined  multiple  shock  interaction,  Waltrup 
and  Billig32  investigated  the  flow  at  Mach  numbers  in  the  range  1.53  <  Mu  <  2.72,  various  unit 
Reynolds  numbers  ranging  from  2.2xl07  <  Re/m  <  6.3xl07,  and  confinement  levels  ranging 
from  0.005  <  9u/r  <  0.020.  At  Mu  =  1.53  they  observed  a  single  normal  shock.  At  all  other 
Mach  numbers  tested,  they  reported  finding  an  oblique  shock  system  and  postulated  that  the 
shock  system  was  a  series  of  weak  reflected  oblique  shocks  with  totally  supersonic  flow 
following  each  shock. 

The  Reynolds  number  effect  on  the  shock  structure  is  clearly  seen  by  comparing  the 


results  of  Mateer,  et  al.33  (Re/m  =  2.3xl07,  Mu  =  1.5,  Sy/r  =  0.202),  who  found  a  single 
normal  shock,  to  those  of  Om  and  Childs35  (Re/m  =  4.9x1 06,  Mu  =  1.49,  Sy/r  =  0.198)  who 
found  a  repeated  normal  shock  system.  Increasing  the  Reynolds  number  for  almost  identical 
Mach  numbers  and  confinement  parameters  causes  the  shock  structure  to  shift  from  a  multiple 
to  a  single  shock  system.  The  effect  of  Mach  number  on  the  shock  structure  is  not  clear,  but  a 
tendency  toward  repeated  shock  systems  as  Mach  number  increases  is  evidenced  by  the 
increase  in  the  length  of  the  pressure  rise  with  increasing  Mach  number.  The  results  of 
Waltrup  and  Billig32  suggest  that  a  single  normal  shock  may  not  exist  in  high  Mach  number 
confined  flows,  even  with  small  However,  the  results  of  Cuffel  and  Back36  indicate  that 
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a  single  normal  shock  can  exist  in  a  high  Mach  number  flow  with  relatively  large  5u/r.  A 
comparison  of  these  two  works  also  contradicts  the  expected  Reynolds  number  effect 
( increasing  Re/m  leading  to  a  single  shock  system).  The  cause  of  this  apparent  discrepancy  is 
not  clear,  indicating  the  need  for  further  work  to  establish  what  types  of  shock  structures  can 
be  expected  at  higher  Mach  numbers,  i.  e.  Mach  numbers  of  approximately  1.8  and  above. 

In  reviewing  the  body  of  experimental  work  related  to  the  internal  shock  wave/turbulent 
boundary  layer  interaction,  several  areas  requiring  further  study  are  evident.  Further  work  is 
needed  to  determine  when  the  transition  from  a  single  normal  shock  to  multiple  normal  shocks 
to  repeated  oblique  shocks  occurs  at  various  Mach  numbers.  Reynolds  numbers,  and 
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confinement  parameters.  The  exact  nature  of  the  repeated  oblique  shock  system  is  also  unclear. 
From  the  existing  data  it  is  difficult  to  determine  if  the  oblique  shocks  are  strong  (subsonic 
flow  after  each  shock)  or  weak  (supersonic  flow  after  each  shock).  The  steadiness  of  the 
interaction  in  a  constant  area  duct  is  also  questionable.  Small  scale,  high  frequency  motion  of 
the  shock  would  lead  to  false  turbulence  measurements.  This  effect  should  be  investigated 
further.  The  topic  of  unsteady  shock  structures  has  received  more  attention  in  the  planar 
studies.  This  may  indicate  that  the  axisymmetric  geometry  is  inherently  steadier  or  it  may 
simply  indicate  that  shock  motion  is  easier  to  detect  in  planar  geometries.  Conspicuously 
absent  are  any  measurements  of  the  turbulence  structure  in  a  repeated  system  of  either  normal 
or  oblique  shocks.  The  amount  of  turbulence  data  available  for  the  confined  single  normal 
shock  interaction  is  also  severely  limited. 

C.  NUMERICAL  WORK 

Numerous  researchers  have  applied  numerical  techniques  for  solving  the  mass  averaged 
Navier-Stokes  equations  to  external  or  unconfined  shock  wave/boundary  layer  interactions, 
including  MacCormack  and  Baldwin37,  Wilcox3**,  Shang  and  Hankey39,  Deiwert40,  Hung  and 
MacCormack41-43,  and  Settles,  et  al.44.  These  studies  all  demonstrate  the  ability  of  shock 
capturing  finite  difference  techniques  to  predict  the  qualitative  features  of  the  interaction.  With 
the  turbulent  calculations,  close  quantitative  agreement,  especially  for  separated  flows  is 
strongly  dependent  on  the  turbulence  model  employed.  In  a  series  of  papers,  researchers  at 
NASA-Ames  have  studied  a  variety  of  shock  wave/turbulent  boundary'  layer  interaction 
flows.--33"34-45'49  Both  confined  and  unconfined  flows  were  considered.  These  calculations 
employed  the  second-order  accurate,  explicit,  time-splitting,  finite  difference  method  of 
MacCormack50  modified  by  MacCormack's  explicit-implicit  characteristics  scheme.51  Several 
turbulence  models  were  employed  including  a  variety  of  algebraic  eddy-viscositv  models,  the 
one-equation  model  of  Rubesin52,  and  the  two-equation  models  of  Jones  and  Launder55  and 
Wilcox  and  Rubesin.54  All  of  the  calculations  showed  good  qualitative  agreement  with 


experimental  results.  Quantitative  agreement  with  experimental  data  was  more  difficult  to 
obtain,  especially  when  boundary  layer  separation  was  present.  Experimental  measurements  of 
surface  pressure  distributions,  skin  friction,  separation  and  reattachment  locations,  mean 
velocity  profiles,  and  turbulence  kinetic  energy  profiles  were  used  to  evaluate  the  various 
computations.  While  none  of  the  models  was  clearly  superior  in  all  respects,  the  higher  order 
models  were  found  to  be  more  universally  applicable  and  typically  offered  better  accuracy  than 
the  lower  order  models.  The  Wilcox-Rubesin  two  equation  eddy  viscosity  model54  emerged  as 
being  superior  for  general  use  in  shock  wave/boundary  layer  interaction  flows.  However,  a 
deficiency  of  the  higher  order,  eddy  viscosity  models  is  that  they  give  poor  predictions  of 
velocity  profiles  downstream  of  the  single  shock  wave/boundary  layer  interaction,  and  this  may 
limit  the  ability  of  the  higher  order  eddy  viscosity  models  to  predict  multiple  shock 
wave/boundary  layer  interactions. 

In  all  of  the  numerical  investigations  mentioned  above,  only  single  shock  interactions 
have  been  calculated.  To  the  author's  knowledge,  no  calculations  of  multiple  normal  or  strong 
oblique  shock  wave/turbulent  boundary  layer  interactions  have  been  performed.  Several 
reported  calculations  of  supersonic  inlet  flows  have  shown  the  ability  to  calculate  single 
terminal  shocks  and  reflected  weak  oblique  shocks  in  confined  ducts  of  various  geometries. 
Fairly  coarse  grids  were  used  in  these  calculations,  thus  the  details  of  the  flow  are  poorly 
resolved.  The  calculations  of  Knight55'57,  Kumar58'6*),  Drummond  and  Weidner61,  Talcott 
and  Kumar8,  Coakley  and  Hsieh62,  and  Hunter,  et  al.65  have  all  employed  some  form  of 
NlacCormack's  finite-difference  method.  Other  techniques  have  been  employed  as  well.  For 
example,  Paynter  and  Chen64  used  a  zonal  analysis  and  Chaussee  and  Pulliam65  used  the 
implicit  finite-difference  method  of  Beam  and  Warming.66 


IIL  NUMERICAL  TECHNIQUES 


The  techniques  used  for  the  numerical  portion  of  this  investigation  are  presented  in  this 
chapter.  The  computational  strategy  was  to  employ  existing  numerical  techniques,  typical  of 
those  being  used  in  industrial  and  research  applications,  and  to  use  existing  software  when 
possible.  The  NASCRIN  computer  code59,  obtained  from  Dr.  Ajay  Kumar  at  NASA-Langlev 
Research  Center,  has  been  extensively  modified  for  the  current  application.  These 
modifications  included  work  with  the  outflow  and  no-slip  boundary  conditions.  A  reference 
plane  characteristics-type  subsonic  outflow  boundary  condition  was  incorporated  and  the  no¬ 
slip  wall  boundary  conditions  were  modified  to  to  use  a  normal  momentum  equation  in 
calculating  wall  pressure.  The  original  code  contained  the  algebraic  Baldwin-Lomax  turbulence 
model67  which  is  a  widely  used  model  in  applied  research  and  development  activities.  As  w.ns 
mentioned  in  the  literature  review,  the  two-equation  Wilcox-Rubesin  turbulence  model54  has 
been  shown  to  be  effective  for  use  in  calculating  shock  wave/boundary  layer  interaction 
flows54-48  giving  slightly  better  performance  than  the  more  widely  used  Jones-Launder 
model55.  For  these  reasons,  the  Wilcox-Rubesin  model  was  selected  for  use  and  was 
incorporated  into  the  modified  NASCRIN  code.  Wall  function  boundary  conditions  for  use 
with  the  Wilcox-Rubesin  model  have  also  been  incorporated.  The  wall  function  method  helps 
avoid  numerical  problems  due  to  the  extremely  stiff  nature  of  the  turbulence  model  equations 
near  the  wall  and  also  offers  a  considerable  savings  in  computer  time.  Finally,  the  NASCRIN 
code  has  been  modified  with  respect  to  data  input  and  output  and  computational  grid 
generation.  In  the  remainder  of  this  chapter,  the  details  of  these  numerical  techniques  are 
given. 

A.  GOVERNING  EQUATIONS 

The  governing  equations  for  the  mean  flow  are  the  same  as  those  used  bv  Kumar.58-59 
i.  e.  the  conservation  form  of  the  mass  averaged,  time-dependent  Navier-Stokes  equations.  A 

18 


generalized  coordinate  transformation  is  applied  to  take  the  equations  from  the  physical  (x,y) 
plane  to  the  computational  (^,rj)  plane  with  a  unit  grid  spacing  being  used  in  the  mutually 
orthogonal  t,  and  rj  directions.  The  coordinate  transformation  is  applied  such  that  the  equations 
remain  in  strong  conservation  form.  The  resulting  equations  are: 

au  aM  aN  . 

"ST  +  -  +  -  =0  (1) 

^  a^  an 


where  the  vectors  U.  M,  and  N  are 
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puu  ♦  y^O,  -  x^Txy 

PVL1  '  xT|°y 

pH  u  •  pu  -  uax-qx )  -  x^tvOy+qy) +  Txy(vyT]-uxT1) 

pv 

Puv  -  y *ox  +  x^Txy 
pvv  -  y4xxy  +  x4oy 

pH  v  -  pv  -  yHuax+qx)  +  x^(voy+qy)  +  txy(-vy^+ux^) 


The  contravariant  velocities  are  defined  as 

(J  -  u/J  =  V  *V^  =  (y^u  -  x^v)/J  (5) 

V  =  v/J  =  V  *Vr|  =(-y^u  +  x^v)/J  (6) 

where  V  is  the  velocity  vector  in  the  (x,y)  plane.  The  contravariant  velocity  components  U 
and  V  represent  the  velocities  perpendicular  to  lines  of  constant  %  and  T),  respectively.  The 
Jacobian  of  the  transformation  is  given  as 
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and  the  metrics  of  the  transformation  are  defined  as 


N=-V 


y^  =  -V 


^  =  TlyJ 

_  dy 


(8) 


The  subscripts  denote  partial  differentiation,  for  example  s  — .  The  total  stresses,  assumed 

dn 


to  be  comprised  of  a  laminar  and  a  turbulent  contribution,  are  defined  by 

ox  =  oxl  4-  CTxl 

Gy  =  Oyl  +  Oyl 

^xy  =  +  ^xyl 

The  laminar  stresses  are  given  by 
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and  the  turbulent  stresses  are  obtained  using  the  Boussinesq  eddy  viscosity  assumption 
described  below.  The  total  heat  flux  is  also  made  up  of  a  laminar  and  turbulent  contribution, 


qx  -  -  Pr,  +  Prt 


m  Vh 


dx 
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All  of  the  derivatives  with  respect  to  x  and  y  are  actually  calculated  in  the  transformed 
coordinates  \  and  r|.  As  an  example,  the  derivative  of  h  with  respect  to  y  is  calculated  from 
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In  these  equations  p,u,v,  and  p  are  the  density,  mean  velocity  in  the  x-direction,  mean  velocity 
in  the  y-direction,  and  static  pressure.  H  is  the  stagnation  enthalpy  and  h  is  the  static  enthalpy. 
The  laminar  viscosity,  |i/,  is  obtained  from  the  Sutherland  viscosity  law  while  the  turbulent 
viscosity,  (it.  is  calculated  from  the  turbulence  models  described  below.  The  laminar  and 
turbulent  Prandtl  numbers  are  taken  as  Pr;  =  0.72  and  Prt  =  0.90,  respectively  In  the  energy 
equation,  the  dependent  variable  is  written  as  pH-p  and  is  physically  interpreted  as  the  total 
internal  energy  per  unit  volume.  One  additional  relation,  beyond  the  turbulence  model  for  |it< 
is  required  to  close  this  set  of  equations.  This  is  taken  as  the  ideal  gas  law 

p  =  pR  T  (18) 

w  here  R  is  the  gas  constant  and  T  is  the  static  temperature. 

Implicit  in  this  formulation  is  the  Boussinesq  eddy  viscosity  assumption.  The  basic 
assumption  is  that  the  apparent  turbulent  stresses,  Tjj1,  which  arise  from  the  mass  averaging 
operation,  are  related  to  the  mean  rate  of  strain  through  a  turbulent  or  eddy  viscosity,  pt-  In 
this  study,  an  isotropic  eddy  viscosity  is  employed  which  has  the  same  numerical  value  for  all 
elements  of  the  stress  tensor  at  a  point,  i.e.  the  eddy  viscosity  is  independent  of  direction.  A 
mathematical  statement  of  the  eddy  viscosity  assumption  is  written  as 


V=  Pui  uj 


dx;  dx 


V  ' J  V 

where  k  is  the  turbulence  kinetic  energy  defined  as 


2  „  3uk 
+  T5ij(4t^— +Pk) 
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k  =  J  uk  uk 


(20) 


The  primes  denote  fluctuating  quantities  and  the  overbars  indicate  a  mass  average  is  taken. 
Equation  (19)  is  the  constitutive  relation  between  the  turbulent  stresses  and  the  mean  rates  of 
strain  that  is  used  to  calculate  the  turbulent  stress  required  in  Equations  (9-11).  When 
expanded,  Equation  (19)  is  identical  in  form  to  Equations  (12-14)  with  p  replaced  by  pk  and 
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the  /  subscript  replaced  by  t.  The  pk  term  represents  an  apparent  pressure  due  to  turbulent 
fluctuations.  The  eddy  viscosity  assumption  does  not  by  itself  constitute  a  turbulence  model. 
Rather,  it  defines  the  framework  in  which  a  turbulence  model  may  be  developed  and  applied. 
All  that  remains  in  the  turbulence  modeling  procedure  is  to  specify  |i.t  and  k  in  terms  of  known 
quantities. 

By  analogy  to  the  kinetic  theory  for  molecular  viscosity,  the  turbulent  viscosity  is 
modeled  by 

Pt  =  pVtLt  (21) 

where  Vt  and  Lt  are  characteristic  velocity  and  length  scales  of  the  turbulent  flow.  In  general, 
these  characteristic  scales  vary  from  point  to  point  in  the  flow  and  may  depend  on  the  history  of 
the  flow  both  upstream  and  downstream  of  the  location  of  interest.  Thus,  a  fundamental 
distinction  exists  between  the  laminar  viscosity,  which  is  a  local  property  of  the  fluid,  and  the 
turbulent  viscosity,  which  is  a  function  of  the  particular  flow  field.  The  two  turbulence  models 
employed  in  this  investigation,  one  due  to  Baldwin  and  Lomax67  and  the  other  to  Wilcox  and 
Rubesin54,  are  based  on  this  framework. 


1.  BALDWIN -LOMAX  TURBULENCE  MODEL 

The  Baldwin-Lomax  model67  is  a  widely  used  two-layer,  algebraic,  isotropic  eddy 
viscosity  turbulence  model.  The  term  algebraic  refers  to  the  type  of  relations  used  to  calculate 
the  turbulent  length  and  velocity  scales.  This  is  a  two-layer  model  in  that  the  turbulent 
viscosity  is  calculated  using  two  distinct  formulations  in  two  distinct  layers,  the  inner  layer 
close  to  a  solid  wall,  and  the  outer  layer  away  from  the  wall.  A  smooth  transition  between  the 
two  formulations  is  enforced  as  follows 


(M-t)inner  y  -  Ycrossover 
(4t)outer  ycrossover  <  y 


(22) 


Here,  y  is  the  normal  distance  from  the  wall  and  the  crossover  point  is  the  closest  location  to 
the  wall  at  which  the  outer  and  inner  turbulent  viscosities  are  equal.  In  the  inner  region,  the 


Prandtl-Van  Driest  formulation  is  used, 


(Pt)inner  -  P  I 


w  here  the  turbulent  length  scale  is 


Lt  =  KyD  , 


the  turbulent  velocity  scale  is  Lt|  col ,  K  is  von  Karman’s  constant,  and  I  col  is  the  magnitude  of 


the  vorticity 


dy  dx 


D  is  the  Van  Driest  damping  factor 


D  =  f  1  -  exp(-y+/A+)] 


y+  _  y\  PwTw 
Pw 


The  outer  turbulent  viscosity  is  calculated  from 


(Pt)outer  =  0:1  CCp  p  F wake  F^klebfy)-  (27) 

In  this  expression,  Cci  is  the  Clauser  constant,  Ccp  is  an  additional  constant,  and  Feebly)  is  the 
Klebanoff  intermittency  function  given  by 


Flkb<y,=[l+5.5^)6]  . 


The  wt 


'mi  is  calculated  from 


F wake  -  ymax  F max 


where  )  A  is  the  y  location  of  the  peak  in  the  function 


F(y)  =  y  I  col  D 


and  FmaA  =  F(ymaA).  The  turbulent  length  scale  in  the  outer  region  is  proportional  to  ymax  and 


the  turbulent  velocity  scale  is  proportional  to  Fma*.  Thus,  the  distribution  of  vorticity  in  the 
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outer  region  is  used  to  establish  these  outer  scales.  Baldwin  and  Lomax67  noticed  in  their 
original  paper  that  the  function  F(y)  can  develop  multiple  peaks  in  the  vicinity  of  a  shock  wave/ 
boundary  layer  interaction.  If  the  wTong  peak  is  selected,  the  function  F^g  will  be  incorrectly 
evaluated.  Visbal  and  Knight  ^  investigated  this  effect  in  more  detail.  Their  results  indicated 
that  two  main  peaks  develop  in  the  separation  region  of  a  shock  wave/boundary  layer 
interaction.  Selection  of  the  inner  peak  causes  artificially  low  values  of  Fwa^e,  while  selection 
of  the  outer  peak  gave  physically  more  reasonable  results  leading  them  to  recommend  use  of 
the  outer  peak.  They  also  noted  that  near  the  separation  point,  where  the  wall  shear  stress  goes 
to  zero,  the  Van  Driest  damping  function  unrealistically  vanishes  causing  low  values  of  the 
computed  turbulent  viscosity.  The  use  of  the  local  shear  stress,  instead  of  the  wall  shear 
stress,  avoided  this  problem.  The  constants  CCp  and  Ckieb  exhibited  a  Mach  number 
dependence.  A  value  of  Ccp=2.08  was  recommended  for  a  Mach  3.0.  nearly  adiabatic  flat  plate 
boundary  layer.  No  specific  recommendation  concerning  C^b  in  a  compressible  flow  was 
made  .  The  original  values  of  the  model  constants  are67 


A+ 
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C,-p  = 

1.6 

Oklcb  = 

0.3 

K 

0.4 

Cd  = 

0.0168 

(30) 


The  kinetic  energy  term,  pk,  in  the  constitutive  relation.  Equation  (19),  is  set  to  zero. 


2.  WILCOX-RUBESIN  TURBULENCE  MODEL 

The  Wilcox-Rubesin  turbulence  model54  is  a  two-equation,  eddy  viscosity  model. 
Two  partial  differential  equations  are  solved  for  the  turbulence  kinetic  energy  and  the  specific 
turbulence  dissipation.  These  two  quantities  are  then  used  in  a  set  of  algebraic  relations  to 
obtain  the  turbulent  length  and  velocity  scales  which  are  combined  to  obtain  a  turbulent 
viscosity.  Wilcox  and  Rubesin54  preferred  to  interpret  the  quantity  k  which  will  appear  in  the 
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turbulence  model  equations  as  a  turbulence  mixing  energy  rather  than  the  total  turbulence 
kinetic  energy.  The  turbulence  mixing  energy  is  related  to  the  fluctuating  velocity  component 
w  hich  is  perpendicular  to  the  mean  flow  direction,  but  was  not  precisely  defined  by  Wilcox  and 
Rubesin^4.  In  this  work,  the  quantity  k  is  taken  as  the  total  turbulence  kinetic  energy.  This 
choice  is  in  agreement  with  that  of  Viegas  and  Horstman4^  jn  their  implementation  of  the 
Wilcox-Rubesin  model.  The  specific  turbulence  dissipation,  to,  defined  as  the  rate  of 
dissipation  of  turbulence  kinetic  energy  per  unit  of  turbulence  kinetic  energy,  is  related  to  the 
dissipation  of  turbulence  kinetic  energy  by 


co  = 


(5*k 


on 


w  here  (3*  is  a  constant  and  e  is  the  dissipation  of  turbulence  kinetic  energy  used  in  the  Jones- 
Launder  model5-\  The  turbulent  viscosity  is  calculated  from 


*  k 

Pt  =  PY  ~ 
co 


(32) 


w  here  k  and  co  are  obtained  from  the  turbulence  model  equations  for  pk  and  pco2  : 


5(pk)  5(pu,k)  ,3(uj)  ,  a  r  .  d(k) 
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dx 


J  J 


=  0 


(33) 


a< pco2)  atpu.co2)  co2  ,  a(u,) 

- +  — 1 - +y-_  t,,1 - 


at 


3xi 


k  ,J  ax, 


(3+2d 


raL,\ 


axk  J  j 


pcoJ 


ax. 


(p/+apt) 


dtio2) 


dx 


J  J 
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(34) 


The  turbulent  lentzth  scale  is 


Li  =  “  . 

CO 


(33) 
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Five  terms  appear  in  each  of  the  turbulence  model  equations  for  k  and  co2,  Equations  (33)  and 
(34).  From  left  to  nght.  these  terms  are  the  time  rate  of  change  followed  by  the  convection, 
production,  viscous  dissipation,  and  viscous  diffusion  terms.  A  viscous  damping  term,  y*.  is 
included  in  the  expression  for  the  turbulent  viscosity  and  is  given  by 

y*  =  [l  -  (1-X2)  exp(-Rei/Rk)]  (36) 

To  allow  a  different  level  of  damping  in  the  co2  equation,  a  second  damping  term,  y,  is  defined 
by 

YY*  =  Yoo[l  -  (1-X2)  exp(-Ret/Roi)]  (36) 

These  viscous  damping  terms  were  required  in  the  model  to  correctly  simulate  the  near  wall 
region  where  molecular  viscosity  dominates.  The  Reynolds  number  of  turbulence  is  expressed 


as 


Ret  =  ^— Ti 


(37) 


The  model  closure  coefficients  are 


10 

Y°°  —  g 


*  =  IT 


Rk  =  l 


R(o  =  2.0 

Viegas  and  Horstman49  found  that  calculations  with  the  original  value  of  y»  =  10/9 
underpredicted  the  wall  shear  stress  in  a  flat  plate  equilibrium  boundary  layer  by  27%  and 
recommended  that  a  value  of  Yoo  =0.90  be  used  instead.  This  set  of  model  equations,  when 
solved  in  conjunction  with  the  mean  flow  equations,  gives  the  turbulent  viscosity,  |it,  and  the 
turbulent  kinetic  energy,  k,  required  in  Equation  (19).  In  Wilcox  and  Rubesin’s54  original 
work,  terms  were  added  to  the  constitutive  relation,  Equation  (19),  to  allow  for  non-aligned 
stress  and  strain.  These  terms,  intended  as  a  correction  to  the  isotropic  model  for  non-isotropic 
effects,  were  found  to  cause  numerical  difficulties  at  times,49  especially  in  calculations  of 
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shock.  wave/boundary  layer  interaction  flows,  and  have  therefore  not  been  included  in  the 
current  formulation. 

B.  NUMERICAL  SOLUTION  METHOD 

The  set  of  governing  equations  given  above  are  solved  using  the  original  explicit,  time- 
dependent,  second-order  accurate,  predictor-corrector,  finite-difference  method  of 
MacCormack69.  A  shock  capturing  approach  on  a  fixed  grid  is  taken.  This  method  has  the 
advantages  of  being  easily  understood  and  implemented  and  is  highly  vectorizable.  In  this 
computational  approach,  the  solution  is  advanced  in  time  from  some  initial  value  plane  to  a  final 
converged,  steady  state  solution.  The  advantage  of  using  a  time-dependent  method  to  obtain 
the  steady  state  solution  is  seen  by  considering  the  nature  of  the  governing  partial  differential 
equations.  The  governing  equations  are  hyperbolic  in  time,  but  display  a  mixed  elliptic- 
hyperbolic  nature  in  the  spatial  coordinates  due  to  the  mixed  subsonic-supersonic  velocities  of 
the  strongly  shocked  boundary  layer  flow  under  consideration.  A  steady  state,  spatial 
marching  technique  would  require  different  computational  procedures  in  the  subsonic  and 
supersonic  portions  of  the  flow.  This  is  avoided  by  using  the  time-dependent  approach  to 
march  the  solution  in  time  until  the  desired  steady  solution  is  obtained,  allowing  the  spatial 
differences  to  be  performed  in  the  same  manner  regardless  of  the  local  flow  velocity. 

In  MacCormack's  method69  the  solution  is  progressed  from  one  time  plane  to  the  next 
in  two  steps,  the  predictor  and  corrector.  The  same  finite-difference  method  is  used  to  solve 
both  the  mean  flow  and  turbulence  partial  differential  equations.  However,  to  illustrate  the 
method  consider  only  the  mean  flow  equations,  Equation  (1),  for  which  the  predictor-corrector 
sequence  is 


Upn+1  =  Un  -  At/A^(D+Mn)  -  At/AqfD+N") 

Ucn+I  =  Un  -  At/A4(D.Mpn+1)  -  At/Ar|(D.Npn+1) 
Un+1  =  l/2(Upn+1  +Ucn+1) 
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The  predictor  step  is  given  by  Equation  (39)  and  the  corrector  step  is  given  by  Equations  (40). 
In  these  expressions,  the  superscript  n  denotes  the  previous  time  plane,  and  the  superscript  n+ 1 
denotes  the  new  time  plane.  The  subscripts  p  and  c  indicate  that  predicted  and  corrected  values 
ot'  the  flow  variables  are  used  to  evaluate  the  vector.  The  time  step  is  At  and  the  gnd  spacings 
are  Ac  and  Arp  The  finite  difference  operators  D+  and  D.  need  further  explanation.  The 
operator  D+.  indicates  that  first  order  accurate  forward  differences  are  used  to  evaluate  the  outer 
spatial  derivatives  in  the  governing  mean  and  turbulence  equations  while  first  order  backward 
differences  are  used  for  all  derivatives  in  the  shear  stress  and  heat  flux  terms.  Equations  ( 12- 
lb.  19 1,  and  for  the  inner  derivatives  of  the  diffusion  terms  in  the  turbulent  kinetic  energy  and 
turbulent  dissipation  rate  equations.  Equations  (33)  and  (34).  The  operator  D_  indicates  the 
opposite  differencing  scheme,  first  order  backward  differences  on  all  outer  derivatives  and  first 
order  forw ard  differences  on  all  inner  derivatives.  In  the  original  MacCormack  scheme,  central 
differences  were  used  on  all  inner  cross-derivatives  in  the  shear  stress  and  heat  flux  terms. 
Central  differences  were  not  used  on  these  terms  in  the  present  code.  This  allows  a  simpler 
algorithm  which  saves  both  computer  memory  space  and  execution  time.  The  use  of  one-sided 
differences  on  these  cross  derivatives  still  retains  second-order  accuracy  for  the  predictor- 
corrector  sequence  and  is  expected  to  have  only  slight,  if  any,  effect  on  the  computed  results.70 
The  sequence  of  operators,  D+,  D.  for  the  predictor-corrector  steps  is  reversed  at  every  time 
step  so  that  the  repeated  sequence  of  operators  is  D+,  D_.  D.,  D+. 

A  fourth-order  numerical  damping  option^9  is  included  in  the  code.  It  is  available  for 
use  when  strong  shocks  are  present  in  the  flow  and  is  required  to  avoid  excessively  large 
numerical  oscillations  in  the  computed  results  near  such  shocks. 

(  .  BOUNDARY  AND  INITIAL  CONDITIONS 

The  numerical  procedure  described  above  is  used  to  solve  the  governing  differential 
equations  at  all  grid  locations  except  the  boundaries.  At  each  new  time  plane  in  the  solution 
procedure,  appropriate  boundary  conditions  must  be  applied.  The  philosophy  taken  in  this 
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work  is  that  while  the  boundary  points  represent  only  a  fraction  of  the  computational  gnd 
points,  the  boundary  conditions  applied  at  these  points  exert  a  large  influence  on  the  computed 
results.  Consequently,  when  formulating  and  applying  the  boundary  conditions,  care  was 
taken  to  realistically  model  the  actual  physical  constraints.  In  the  following  sections,  the 
inflow,  outflow,  solid  wall,  and  symmetry  plane  boundary  conditions  are  treated  separately. 
An  initial  value  plane  must  also  be  specified  to  start  the  calculation.  The  initial  value  plane 
should  be  selected  such  that  large  numerical  transients  which  might  cause  the  solution  to 
become  unstable  are  avoided.  The  choice  of  initial  value  plane  could  also  affect  the  overall  time 
required  for  convergence  to  a  steady  solution  and  could  conceivably  determine  which  steady 
state  solution  is  reached,  since  the  uniqueness  of  the  solution  is  not  assured.  The  actual  initial 
value  plane  used  for  each  calculation  is  stated  where  that  calculation  is  discussed. 

1  INFLOW  AND  SYMMETRY  BOUNDARY  CONDITIONS 

The  inflow  and  symmetry  boundary  conditions  are  treated  together  since  they  are  both 
relatively  simple.  At  the  inflow  boundary  the  flow  is  assumed  to  be  supersonic  except  in  the 
near  wall  regions  where  the  flow  may  be  subsonic.  For  these  conditions  the  governing 
equations  are  either  hyperbolic  or  parabolic  and  one  would  not  expect  information  to  propagate 
upstream.  Accordingly,  the  primitive  variables  u,  v,  T,  p,  k,  and  to  are  held  fixed  at  specified 
incoming  values.  At  a  symmetry  boundary,  the  symmetric  quantities  u,  T,  and  p  are  obtained 
from  a  one  point  (zeroth  order)  extrapolation  (equivalent  to  zero  normal  gradient)  while  the 
antisymmetric  quantity  v  is  set  to  zero.  The  turbulence  quantities  k  and  co  are  also  obtained 
from  one  point  extrapolations. 

2.  OUTFLOW  BOUNDARY  CONDITIONS 

Several  options  are  allowed  at  the  outflow  boundary.  The  turbulence  variables  k  and  co 
are  extrapolated  at  both  subsonic  and  supersonic  outflow  points  using  a  zeroth  order 
extrapolation.  With  supersonic  outflow,  information  is  physically  expected  to  be  swept  out  of 
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the  grid  so  that  either  zeroth  or  first  order  extrapolation  may  be  used  for  u,  v,  p,  and  T.  This 
type  of  extrapolation  is  also  applied  at  the  near  wall  subsonic  points  in  a  parabolic  boundary 
layer  type  flow.  The  reference  plane  characteristic  scheme,  described  below,  may  also  be  used 
for  supersonic  outflow  conditions.  When  the  flow  is  subsonic  at  the  outflow  boundary, 
information  can  propagate  upstream  and  extrapolation  of  the  mean  flow  variables  is  not 
physically  consistent.  As  will  be  shown,  one  of  the  mean  flow  variables  u,  v,  p,  or  T  must  be 
specified  at  the  exit.  Rudy  and  Strikwerda71  show  that  a  proper  treatment  of  a  subsonic 
boundary  is  needed  to  prevent  the  numerical  reflection  of  information  at  the  boundary.  A 
simple  treatment  such  as  extrapolating  u,  v,  and  T  and  holding  p  fixed  can  cause  the  reflection 
of  pressure  waves  back  into  the  computational  grid.  Low  frequency  pressure  waves  can 
actually  be  trapped  in  the  computational  grid.  To  avoid  such  problems,  a  non-reflective, 
reference  plane  characteristics  scheme,  similar  to  that  of  Cline72*73  has  been  implemented. 

A  method  of  characteristics  analysis  of  the  non-conservative  form  of  the  governing 
equations  is  performed  to  find  the  characteristic  and  compatibility  equations.  The  £  direction  is 
perpendicular  to  the  boundary  and  the  T|  direction  is  parallel  to  the  boundary.  All  r\  derivatives 
and  viscous  terms  are  evaluated  using  finite  differences  and  are  treated  as  source  terms.  The 
resulting  set  of  compatibility  and  associated  characteristic  relations  are 

dp  -  Qjdp  =  H'adt  along  ^=C/  (41) 

f^2du  +  Q3dv  =  (Q2^2  +  ^3^3)dt  along  f/  (42) 

Q4du  -  ftsdv  +  dp  =  (fii'P i  +  SVP2  -  ^5^3  +  ^4)dt 


and 

-Q^du  .  Q-jdv  +  dp  =  (fli^i  -  0^2  +  ^7^3  +  'Tridt 


The  Q,  and  VP,  are  defined  in  Appendix  A.  The  contravariant  velocity  U  was  previously 
defined  in  Equation  (5)  and  a  is  the  speed  of  sound  defined  as 
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a  =  V  yRT  . 


(45) 


The  total  denvatives  in  Equations  (41-44)  are  applied  along  the  characteristics.  The  three 
characteristic  directions  are  obtained  from  the  equations  for  d^/dt.  When  the  flow  is  supersonic 
in  the  q,  direction,  all  three  characteristics  extend  backwards  into  the  computational  domain 
from  the  boundary,  indicating  that  the  values  of  u,  v,  p,  and  T  at  the  exit  are  determined  from 
the  upstream  flow.  When  the  flow  is  subsonic  in  the  q  direction,  the  first  two  characteristics 
extend  into  the  domain  while  one  extends  out  of  the  computational  domain,  such  that  one  of  the 
flow  v  ariables  must  be  specified.  At  an  inflow  boundary,  this  same  analysis  shows  that  four 
flow  variables  must  be  specified  for  supersonic  inflow,  while  only  three  flow  variables  should 
be  specified  for  subsonic  inflow. 

3.  SOLID  WALL  BOUNDARY  CONDITIONS 

Adiabatic,  no  slip,  solid  wall  boundary  conditions  are  enforced.  The  method  of 
applying  these  boundary  conditions  depends  on  the  turbulence  model  employed.  With  the 
Baldwin-Lomax  model  the  mean  flow  equations  are  integrated  to  the  wall  with  the  first  grid 
point  off  the  wall  placed  in  the  viscous  sublayer  to  provide  sufficient  resolution  of  the 
extremely  steep  velocity  gradients  in  the  near  wall  region.  The  u  and  v  velocity  components  Eire 
set  to  zero  at  the  wall.  The  wall  temperature  is  obtained  from  a  zeroth  order  extrapolation 
which  is  equivalent  to  enforcing  a  zero  normal  temperature  gradient.  The  wall  pressure  is 
obtained  from  a  normal  momentum  equation. 

Two  options  are  available  for  applying  the  wall  boundary  conditions  when  the  Wilcox- 
Rubesin  model  is  used.  The  first  is  to  integrate  the  governing  equations  to  the  wall  as  is  done 
with  the  Baldwin-Lomax  model.  However,  a  much  finer  grid  spacing  is  required  to  resolve  the 
gradients  in  the  turbulence  quantities.  With  the  Baldwin-Lomax  model  a  value  of  v*  between  1 
and  5  at  the  first  grid  point  off  the  wall  is  sufficient.  To  accurately  solve  the  Wilcox-Rubesin 
turbulence  model  equations,  this  value  of  y+  must  be  reduced  by  a  factor  of  ten  with  a  value 
below  0. 1  recommended.  The  mean  flow  variables  u,  v,  p,  and  T  are  obtained  as  for  the 
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Baldwin-Lomax  model.  The  turbulent  kinetic  energy  is  set  to  zero  at  the  wall.  Physically,  the 

specific  turbulence  dissipation  goes  to  infinity  at  the  wall.  However,  the  following  asymptotic 

value  of  to  is  recommended  by  Wilcox  and  Rubesin54  for  a  smooth  impervious  w  all 
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and  was  found  by  Viegas  and  Horstman4^  to  be  valid  up  to  y+=10. 

Initial  computations  using  Equation  (45)  and  integration  to  the  wall  showed  that  the 
turbulence  model  equations  were  extremely  stiff  in  the  near  wall  region  causing  numerical 
problems  for  even  a  simple  flat  plate  boundary  layer  case.  Accurate  solutions  could  be 
obtained  but  required  inordinate  amounts  of  computer  time  due  to  the  extremely  fine  gnd 


spacing. 

In  an  effort  to  avoid  these  problems,  Viegas  and  Rubesin74  develop'd  a  wall  function 
boundary  condition  for  use  with  the  Jones-Launder  and  Wilcox-Rubesin  turbulence  models 
that  is  applicable  to  adiabatic,  no-slip  wall  conditions.  Their  results  were  most  promising  with 
the  general  conclusion  that  the  wall  functions  improved  the  agreement  of  computed  results  with 
experiments,  especially  in  the  regions  of  shock  wave/boundary  layer  interactions  and  also 
helped  avoid  the  stiffness  problems  associated  with  the  near  wall  region.  This  formulation 
required  that  the  first  grid  point  off  the  wall  be  located  in  the  fully  turbulent  log  region  with  a 
y+  of  at  least  40.  The  upper  limit  on  y+  at  the  second  grid  point  off  the  wall  was  constrained  to 
be  less  than  5+/7.  Thus,  a  vast  reduction  in  computer  time  could  be  achieved.  In  a  later  paper. 
Viegas,  et  al.75  presented  an  extended  version  of  the  wall  function  method  with  improved 
performance  in  the  calculation  of  separated  flows.  The  new  method  was  also  applicable  to 
non-adiabatic  walls  and  included  a  formulation  for  use  in  the  viscous  sublayer.  The  original 
wall  function  method  was  incorporated  into  the  NASCRIN  code  and  was  used  for  a 
preliminary  set  of  calculations  with  adiabatic  wall  conditions.  For  a  flat  plate  boundary  layer 
case,  the  wall  temperature  was  found  to  increase  with  downstream  distance,  eventually 
exceeding  the  incoming  stagnation  temperature.  The  reason  for  this  was  not  apparent. 
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Although  the  wall  function  method  applied  in  the  log  layer  appeared  to  be  most 
promising,  the  problem  with  the  wall  temperature  was  disturbing.  For  the  current  applications 
where  detailed  experimental  measurements  through  the  boundary  layer  are  going  to  be 
evaluated  and  compared  to  numerical  predictions,  a  finer  grid  spacing  than  the  log  layer  wall 
function  method  provides  is  desired.  For  this  reason,  the  wall  function  formulation  for  use  in 
the  viscous  sublayer75  is  used  in  this  study.  With  this  formulation,  the  grid  spacing 
requirements  near  the  wall  are  the  same  as  with  the  Baldwin-Lomax  model.  Hence,  the  two 
techniques  may  be  compared  on  the  same  grid.  The  sublayer  formulation  is  also  much  easier  to 
vectorize  than  the  log  layer  method  while  it  still  offers  the  advantage  of  reducing  the  numerical 
stiffness  problems.  The  sublayer  type  wall  function  approach  is  outlined  below. 

At  the  wall  u  =  v  =  qy  =  0  and  the  vector  N  at  the  wall  becomes 


0 

-y^Pw  + 

■y^^xy^w  +  x£Pw 


(46) 


It  is  this  vector  that  is  needed  in  Equation  (1)  to  impose  the  wall  boundary  conditions.  The 
wall  shear  stress,  (ixy)w.  the  wall  pressure,  pw,  and  the  wall  temperature  Tw  are  required  to 
evaluate  this  expression.  The  wall  pressure  is  obtained  from  the  same  normal  momentum 
equation  used  with  the  Baldwin-Lomax  model.  In  the  following  development  of  the  wall 
function  relationships,  the  coordinate  system  is  slightly  different  than  in  the  mean  flow 
equations.  The  y  direction  is  perpendicular  to  the  wall  and  positive  away  from  the  wall.  The  u 
velocity  is  the  velocity  along  the  wall  and  the  heat  flux  and  shear  stress  terms  are  defined  by 
gradients  perpendicular  to  the  wall.  The  first  grid  point  off  the  surface  is  denoted  by  a 
subscript  2  and  the  edge  of  the  viscous  sublayer  is  denoted  by  a  subscript  v. 

The  wall  temperature  is  obtained  from  the  reduced  energy  equation  in  the  viscous 
sublayer 

q  =  qw  +  ux 


it 


(47) 


where  q  is  the  heat  flux  normal  to  the  wall,  u  is  the  velocity  parallel  to  the  wall  and  t  is  the 
shear  stress  aligned  with  wall.  For  an  adiabatic  wall,  equation  (47)  may  be  written  as 
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which  is  integrated  to  yield 


(49) 
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This  expression  is  applied  at  point  2,  where  u  =  U2,  T  =  T2  to  give  the  wall  temperature,  Tw. 

The  wall  shear  stress  is  obtained  from  the  reduced  momentum  equation  in  the  sublayer, 

3u  dp 

i— =  Tw  + 
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The  viscosity  is  assumed  to  vary  linearly  with  temperature,  p  =  pw  (T/Tw).  Using  this 
assumption  and  the  temperature  distribution  given  by  Equation  (49),  Equation  (50)  is  integrated 
to  give  an  expression  for  the  wall  shear  stress 
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(51) 
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Equations  (49)  and  (51)  are  used  to  evaluate  the  wall  shear  stress  and  temperature  given  values 
of  u,  T,  and  y  at  point  2  and  dp/dx  at  the  wall,  the  only  restriction  being  that  this  grid  point  be 
located  in  the  viscous  sublayer,  yj  <  yv.  The  edge  of  the  viscous  sublayer  is  taken  as  the 
location  where 

yvV  Pvky/pw  (5?) 
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A  parabolic  kinetic  energy  distribution 

k  =  k, 
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is  assumed  in  the  viscous  sublayer.  Equations  (52)  and  (53)  are  combined  to  give  an 
expression  for  the  location  of  the  viscous  edge 
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Special  conditions  must  also  be  applied  to  the  turbulence  model  equations  when  employing  the 
wall  function  method.  From  the  kinetic  energy  distribution,  equation  (53),  one  sees  that  a  zero 


normal  gradient  in  k  exists  at  the  wall,  i.  e. 


St  =0 


This  condition  is  imposed  on  the  normal  derivatives  of  k  at  the  wall  which  appear  in  the 

turbulence  model  equation  for  k.  The  turbulence  kinetic  energy  production  and  dissipation 

terms  also  receive  special  treatment.  The  production  of  turbulence  kinetic  energy  at  point  2  is 

set  to  zero.  The  turbulence  dissipation  is  set  using  Equation  (45), 

20Pw 

“2  -  r-^2  (56) 
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and  the  dissipation  term  in  the  kinetic  energy  equation  is  then  calculated  using  this  value  of  o>2. 
By  specifying  con,  the  need  to  solve  Equation  (34)  at  point  2  is  eliminated. 

A  brief  summary  of  the  wall  wall  function  boundary  condition  is  now  given.  Equations 
(49)  and  (51)  are  solved  for  the  wall  temperature  and  wall  shear  stress  needed  in  the  vector  Nw. 
The  wall  pressure  is  obtained  from  a  normal  momentum  equation.  The  turbulence  dissipation 
is  set  using  Equation  (56)  and  the  turbulence  kinetic  energy  is  set  to  zero.  Additionally,  the 
condition  of  zero  gradient  in  the  turbulence  kinetic  energy  as  stated  in  Equation  (55)  is  imposed 
when  evaluating  the  turbulence  model  equations  at  point  2.  These  relations  are  used  at  solid 
wall  boundary  points  and  have  been  found  to  allow  a  coarser  grid  spacing  near  the  wall  which 
helps  avoid  numerical  problems  encountered  when  the  turbulence  equations  are  integrated  to 
the  wall  and  allows  a  larger  time  step. 


D.  TEST  CALCULATIONS 

Several  test  calculations  were  performed  to  validate  the  computational  procedure 
outlined  in  the  preceding  sections.  All  of  the  computations  reported  here  were  performed  on 
the  NCSA  CRAY-XMP/48  vector  processing  supercomputer.  This  computer  is  configured 
such  that  each  of  the  four  available  processors  act  independently,  i.e.  no  parallel  processing, 
and  the  core  memory  of  this  computer  was  sufficiently  large  to  execute  all  of  the  runs  in  the 
core  memory  without  having  to  resort  to  the  use  of  the  Solid  State  Disk  (SSD)  device  which  is 
available.  The  SSD  is  essentially  a  very  large  and  fast  RAM  disk  device.  Three  test  cases  were 
executed.  The  first  was  a  Mach  1.48  equilibrium  flat  plate  turbulent  boundary  layer,  the 
second  was  a  Mach  2.94  equilibrium  flat  plate  turbulent  boundary  layer,  and  the  third  was  a 
Mach  2.94  oblique  shock/turbulent  boundary  layer  interaction  at  a  20  deg.  compression  corner. 
The  results  of  these  three  cases  are  now  considered  separately. 

1.  TEST  CASE  1 

This  test  case  was  performed  mainly  as  a  check  for  the  wall  function  boundarx 
conditions  but  is  also  a  useful  check  on  the  proper  performance  of  the  main  computational 
scheme  and  the  Wilcox-Rubesin  turbulence  model.  Computations  were  made  of  an  adiabatic 
equilibrium  flat  plate  boundary  layer  at  M  =  1.48,  To  =  283  K,  and  Re§  =  5.83  x  105,  5  =  25 
mm.  These  flow  conditions  correspond  to  the  experimental  conditions  of  Mateer,  et  al.7\ 
These  experiments  were  performed  in  a  constant  area  duct  of  radius  0.1239  m  (4.875  in.). 
This  case  was  also  computed  by  Viegas  and  Rubesin74  using  both  integration  to  the  wall  and 
log-layer  type  w'all  function  boundary  conditions.  The  integration  to  the  w'all  employed  the 
Wilcox-Rubesin  model  while  the  reported  results  for  the  wall  function  calculations  employed 
the  Jones-Launder  model. 

To  compute  this  case  the  computational  grid  was  set  up  with  a  symmetry  plane  located 
0.12m  from  the  flat  plate.  In  the  direction  normal  to  the  wall,  52  grid  points  were  used.  A 
geometric  progression,  with  a  grading  ratio  of  1.05,  was  used  to  stretch  the  grid  near  Te  wall. 


The  geometric  progression  was  applied  to  the  first  41  points  off  the  wall  and  a  constant  grid 
spacing  was  used  for  the  remaining  points.  The  spacing  between  the  first  point  and  the  wall 
was  \'2  =  6.0  x  10-4  m  (vs*  =  170).  resulting  in  21  grid  points  being  located  in  the  boundary 
layer.  In  the  flow  direction  a  constant  grid  spacing  of  Ax/5  =  1.0  was  used.  The  incoming 
boundary  layer  for  the  calculation  was  obtained  using  the  Baldwin-Lomax  turbulence  model  to 
calculate  a  boundary  layer  growing  from  a  uniform  initial  profile.  The  results  from  the 
Baldwin-Lomax  calculation  were  used  to  start  the  Wilcox-Rubesin  calculation  by  assuming  a 
uniform  turbulence  intensity  through  the  boundary  layer  to  give  the  turbulence  kinetic  energy, 
k.  The  specific  dissipation,  co,  was  obtained  using  this  value  of  k,  the  turbulent  viscosity  and 
density  from  the  Baldwin-Lomax  calculation,  and  Equation  (31)  with  y*  =  1.0.  The  calculation 
was  then  continued  using  the  log-layer  wall  functions  and  the  Wilcox-Rubesin  model. 
Inaccuracies  in  the  profiles  of  the  turbulence  quantities  caused  by  this  crude  patching  together 
of  the  two  models  were  found  to  quickly  damp  out  a  few  grid  locations  downstream.  The 
initial  value  plane  was  obtained  by  applying  the  incoming  flow  variable  profiles  to  all 
downstream  locations.  Extrapolation  was  used  at  the  outflow  boundary  and  symmetry 
conditions  were  imposed  at  the  freestream  boundary.  The  original  constants  (Equation  38)  for 
the  Wilcox-Rubesin  model  were  used  with  the  exception  of  y«,=  0.90  being  used,  as  this  value 
of  yx  was  the  one  used  by  Viegas  and  Rubesin74.  Numerical  damping  was  not  used.  The 
calculation  w  as  made  in  segments  by  dividing  the  flow  direction  into  three  pieces.  The  final 
calculation  was  made  on  a  52  x  51  grid  and  was  nin  for  12000  iterations  requiring  323  cpu 
seconds  giving  a  computational  rate  of  1.02  x  1(L5  seconds  /  grid  point  /  time  step.  Thi  w  as 
sufficient  physical  time  to  allow  a  panicle  in  the  freestream  to  traverse  the  computational 
domain  5  times. 

To  evaluate  the  computed  results  a  companson  to  Viegas  and  Rubesin's74  integration  to 
the  wall  is  made  in  Figures  5  and  6.  In  Figure  5  the  computed  u  velocity  profiles  are  shown 
with  good  agreement  between  the  two  computations.  A  more  severe  test  is  shown  in  Figure  6 
where  the  computed  turbulent  kinetic  energy  profiles  are  compared.  Again,  good  agreement  is 


seen.  These  results  were  most  encouraging  as  they  indicate  that  the  current  implementation  of 
the  Wilcox-Rubesin  model  with  the  log-layer  wall  function  boundary  conditions  could 
accurately  reproduce  calculations  made  with  an  entirely  independent  code.  However,  upon 
further  study  of  the  computed  results,  the  wall  temperature  was  not  found  to  be  constant.  The 
wall  temperature  was  seen  to  increase  with  downstream  distance  eventually  exceeding  the 
incoming  stagnation  temperature.  This  condition  is  physically  unrealistic  for  this  adiabatic  wall 
calculation.  Similar  calculations  employing  the  Wilcox-Rubesin  model  and  the  sublayer  wall 
functions  did  not  display  this  undesirable  trend.  The  reason  for  this  discrepancy  is  unknown 
and  because  of  this  only  the  sublayer  type  wall  functions  are  employed  in  the  remainder  of  this 
work. 

2.  TEST  CASE  2 

For  the  second  test  case  the  Mach  2.94  adiabatic  flat  plate  equilibrium  boundary  layer 
investigated  experimentally  by  Kuntz84  was  computed  using  both  the  Baldwin-Lomax  and 
Wilcox-Rubesin  turbulence  models.  The  velocity  measurements  of  Kuntz84  were  made  using 
the  same  two-component,  coincident  LDV  system  to  be  used  in  this  study.  These 
measurements  were  made  on  the  floor  of  a  10.2  x  10.2  cm2  wind  tunnel  with  the  upstream 
stagnation  conditions  being  po  =  483  kPa  (70  psia)  and  To  =  303°  K.  The  boundary  layer 
thickness  was  5  =  8.2  mm. 

The  same  computational  grid  in  the  y-direction  was  used  for  both  the  Baldw'in-Lomax 
and  the  Wilcox-Rubesin  calculations.  Fifty-five  grid  points  were  used  in  the  y-direction, 
normal  to  the  wall.  These  grid  points  were  distributed  into  4  regions.  A  geometric  progression 
was  applied  in  the  first  three  regions  away  from  the  wall.  In  the  first  region,  closest  to  the 
wall,  20  points  were  used  and  the  grading  ratio  was  1.15.  The  location  of  the  first  grid  point 
away  from  the  wall  was  at  y2  =  1.0  x  10'^  (y2+  =  2.5).  In  the  next  region,  20  more  points 
were  used  with  a  grading  ratio  of  1 .05.  In  the  third  region,  10  points  wrere  used  with  a  grading 
ratio  of  1.30.  Equal  grid  spacing  was  used  for  the  remaining  points  in  the  fourth  region.  The 
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freestream  boundary  was  located  0.1  m  from  the  wall.  The  initial  boundary  layer  profiles  and 
initial  rime  planes  were  obtained  using  the  same  procedure  described  above  for  the  Mach  1.48 
boundary  layer  calculation  and,  as  before,  the  calculation  was  performed  in  segments. 

For  the  final  Baldwin-Lomax  calculation,  41  grid  points  were  used  in  the  flow 
direction  with  a  constant  spacing  of  0.015  m  with  Ax/6  =  1.8  at  the  location  of  interest.  Visbal 
and  Knight's6^  recommended  value  of  Ccp  =  2.08  was  used;  otherwise  the  original  model 
constants  in  Equation  (30)  were  used.  No  slip  boundary  conditions  were  applied  at  the  wall, 
symmetry  boundary  conditions  were  applied  at  the  freestream  boundary,  one  point 
extrapolation  was  used  at  the  exit,  and  all  incoming  flow  variables  were  held  fixed.  The 
explicit-implicit  option  available  in  the  NASCRIN  code  was  used  with  a  Courant  number  of  3. 
No  damping  was  used.  Special  care  was  taken  to  ensure  that  the  absolute  maximum  outer  peak 
in  the  function  F(y)  was  selected.  Several  local  peaks  in  this  function  were  observed.  The 
solution  was  run  for  30,000  iterations  requiring  499  cpu  seconds  giving  a  computational  rate  of 
7.4  x  10-6  seconds  /  grid  point  /  time  step.  A  panicle  in  the  freestream  w'ould  have  had 
sufficient  time  pass  through  the  computational  domain  1.4  times.  Thus,  this  solution  may  not 
be  totally  converged. 

The  final  Wilcox-Rubesin  computation  was  performed  on  a  gnd  with  21  points  in  the 
flow  direction  and  an  equal  spacing  of  0.015  m  giving  Ax/5  =  1.8  at  the  location  of  interest. 
The  patching  between  the  Baldwin-Lomax  and  the  Wilcox-Rubesin  calculations  was  performed 
at  a  location  0.75  m  upstream  of  the  initial  column  on  the  final  gnd.  Sublayer  type  wall 
functions  were  used  at  the  solid  wall,  otherwise  the  boundary  conditions  were  the  same  as  in 
the  Baldwin-Lomax  calculation.  A  Courant  number  of  0.95  was  used  with  the  totally  explicit 
MacCormack  scheme  and  the  model  constants  used  were  those  given  in  Equation  (38)  with  the 
exception  of  .  The  solution  was  run  for  45,000  iterations  using  Yo°  =  9/10,  giving  the  first 
solution.  The  solution  was  then  advanced  for  20.000  more  iterations  with  y°o  =  10/9,  giving 
the  second  solution.  After  the  total  of  65,000  iterations  a  particle  in  the  freestream  would  have 


had  sufficient  time  to  pass  4.25  times  through  the  computational  domain.  The  computational 
rate  was  1.39  x  10'5  sec  /  grid  point  /  time  step. 

A  comparison  with  the  experimental  data  will  now  be  made  for  these  three 
computations.  The  computed  values  of  the  wall  shear  stress  for  the  Baldwin-Lomax,  700=  9/10 
Wilcox-Rubesin,  and  7*,=  10/9  Wilcox-Rubesin  calculations  are  120,  113  and  99  N/m~, 
respectively.  These  correspond  to  friction  velocities  of  uT  =  25,  24.2,  and  22.7  m/s.  The 
experimental  value  of  the  friction  velocity  was  obtained  from  the  velocity  measurements  using  a 
wall-wake  profile  curve  fit  to  the  data  yielding  a  value  of  ut  =  24.6  m/s.  Thus,  all  three 
computations  give  reasonable  values  for  the  wall  shear  stress  with  the  700=  10/9  Wilcox- 
Rubesin  calculation  underpredicting  the  wall  shear  stress  somewhat.  This  trend  is  in  agreement 
with  the  results  of  Viegas  and  Horstman4^.  Streamwise  velocity  profiles  for  the  three 
computations  are  compared  to  the  experimental  results  in  Figure  7.  All  three  computations  do  a 
fair  job  of  predicting  the  velocity  profile.  The  Baldwin-Lomax  calculation  tends  to  overpredict 
the  velocity  over  the  entire  boundary  layer.  Both  Wilcox-Rubesin  calculations  give  better 
predictions  near  the  boundary  layer  edge  but  also  overpredict  the  velocity  closer  to  the  wall. 
The  turbulence  kinetic  energy  profiles  for  the  two  Wilcox-Rubesin  calculations  are  compared 
to  the  data  in  Figure  8.  Both  calculations  overpredict  the  turbulence  kinetic  energy  within  the 
boundary  layer  with  the  700  =  10/9  solution  giving  closer  agreement  to  the  data.  In  the 
freestream,  both  calculations  give  zero  turbulence  kinetic  energy  while  the  experimental  results 
indicate  a  freestream  turbulence  intensity  is  present.  Part  of  this  freestream  turbulence  may  be 
attributed  to  the  limited  resolution  of  the  counter  clocks  used  for  the  LDV  measurements.  In  a 
high  speed  flow,  this  clock  resolution  problem  leads  to  an  overprediction  of  the  freestream 
turbulence  intensity.  The  shape  of  the  two  computed  profiles  are  very  similar  displaying  a  peak 
in  the  turbulence  kinetic  energy  very  close  to  the  wall  and  a  second  peak  in  the  middle  of  the 
boundary  layer.  This  second  peak  is  also  seen  to  occur  in  the  experimental  results. 
Measurements  were  not  made  close  enough  to  the  wall  to  reveal  the  inner  peak  in  the  turbulence 
kinetic  energy.  The  turbulence  kinetic  energy  is  a  very  sensitive  variable  and  the  differences 
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between  the  computed  and  experimental  results  are  not  considered  to  be  excessive.  However, 
the  overprediction  of  the  turbulence  kinetic  energy  could  help  explain  the  slightly  fuller 
computed  mean  velocity'  profiles. 

3.  TEST  CASE  3 

For  the  third  test  case  a  separated  oblique  shock/turbulent  boundary  layer  interaction 
was  selected.  The  experimental  measurements  of  Kuntz84  for  a  Mach  2.94  oblique 
shock/turbulent  boundary  layer  interaction  at  a  20°  compression  comer  are  compared  to 
computed  results.  Calculations  employing  the  Baldwin-Lomax  model  have  been  completed 
and  calculations  with  the  Wilcox-Rubesin  model  are  currently  being  performed.  A  description 
of  the  Baldwin-Lomax  calculation  follows.  The  grid  in  the  y-direction  is  the  same  as  that  used 
above  for  the  Mach  2.94  equilibrium  boundary  layer  computation.  Results  from  the  calculation 
above  were  also  used  to  specify  the  incoming  flow  variables.  The  boundary  conditions  and 
turbulence  model  constants  are  the  same  as  above  but  now  the  totally  explicit  MacCormack 
scheme  is  used  with  a  Courant  number  of  0.95.  The  computational  domain  extended  from  x  = 
0  to  x  =  0.0825  m  with  the  20°  ramp  beginning  at  x  =  0.040  mm.  Two  grids  were  employed, 
varying  only  in  the  streamwise  grid  spacing.  For  the  first  grid,  34  points  were  used  with  a 
constant  spacing  of  Ax  =  2.5  mm,  Ax/5o  =  0.3.  For  the  second  grid  45  points  were  used  with 
the  grid  spacing  varying  from  Ax/5o  =  0.49  upstream  of  the  separation  point,  to  Ax/5o  =  0.12 
in  the  separation  region,  to  Ax/5o  =  0.31  downstream  of  reattachment.  The  first  calculation, 
(a),  employed  the  34  point  grid  and  was  run  for  30,000  iterations  using  a  damping  coefficient 
of  0.05.  This  calculation  was  extended  for  30,000  more  iterations  using  a  damping 
coefficient  of  0.20  to  give  the  second  calculation,  (b).  The  third  calculation,  (c),  was 
performed  on  the  refined  grid  using  a  damping  coefficient  of  0.20  and  was  run  for  16,000 
iterations.  For  computations  (a)  and  (b)  sufficient  physical  time  was  allowed  for  a  particle  in 
the  freestream  to  pass  through  the  domain  5  times  while  in  (c)  a  particle  had  time  to  pass 
approximately  2.5  times  through  the  domain.  However,  in  (c)  the  separation  location  was  not 


found  to  change  significantly  for  1000  iterations.  In  all  three  cases  the  computational  rate  v.as 
approximately  7.8  x  Iff6  sec  /  grid  point  /  time  step. 

Results  for  these  three  calculations  are  shown  in  Figure  8  where  the  wall  static  pressure 
non-dimensionalized  by  the  freestream  static  pressure  has  been  plotted  with  the  downstream 
distance  shifted  such  that  the  comer  occurs  at  an  abscissa  value  of  zero.  All  three  calculations 
predicted  a  small  separation  region  at  the  comer  but  gready  underpredicted  the  upstream  extent 
of  the  separation.  This  result  is  typical  of  the  Baldwin-Lomax  model. Oscillations  in  the 
pressure  are  seen  in  case  (a)  indicating  the  need  for  increased  numerical  damping.  Cases  tbi 
and  (ci  are  nearly  identical  indicating  that  grid  independence  has  been  achieved,  at  least  in  the 
flow  direction.  The  effects  of  the  grid  spacing  normal  to  the  wall  need  to  be  investigated 
further.  LDV  measurements  of  the  velocity  field  are  also  available  for  this  configuration  but 
have  not  yet  been  compared  to  the  computations.  As  was  mentioned  above,  calculations  with 
the  Wilcox-Rubesin  model  are  currently  being  made  for  this  test  case  and  are  expected  to  give 
much  improved  results. 


IV.  EXPERIMENTAL  SETUP 


The  experimental  facilities  and  procedures  used  during  this  investigation  are  described 
in  this  chapter.  In  part  A  attention  is  focused  on  the  air  flow  facility  and  small  scale  wind 
tunnel/test  section  used.  The  preliminary  experiments  performed  to  check  the  operating 
characteristics  of  the  wind  tunnel  are  also  described.  In  pan  B  details  of  the  procedures  and 
equipment  used  to  perform  the  Schlieren,  surface  flow  visualization,  and  the  wall  static 
pressure  measurements  are  given.  Finally,  in  pan  C  the  laser  Doppler  velocimeter  (LDV) 
system  is  discussed.  This  system  will  be  used  in  the  remainder  of  this  study  to  investigate  the 
Mach  1.6  shock  wave/turbulent  boundary  layer  interaction.  LDV  operating  proceedures  and 
data  reduction  and  analysis  considerations  are  given. 

A.  EXPERIMENTAL  AIR  FLOW  FACILITY  AND  TEST  SECTION 

The  experimental  portion  of  this  study  was  performed  in  the  Gas  Dynamics  Laboratory 
of  the  Department  of  Mechanical  and  Industrial  Engineering.  This  high  pressure  air  How 
facility  consists  of  two  compressors  which  supply  air  to  a  storage  tank  farm  with  an 
approximate  volume  of  140  m3.  The  first  compressor,  built  by  Ingersoll-Rand,  delivers  41 
kg/min  at  960  kPa  ( 1200  SCFM  at  125  psig),  and  the  second,  a  Gardner-Denver  compressor, 
supplies  20  kg/min  at  760  kPa  (600  SCFM  at  95  psig).  Thus,  at  pressures  below  760  kPa,  a 
continuous  flowrate  of  61  kg/min  is  possible.  A  six  inch  supply  line  connects  the  tank  farm  to 
a  large  axisymmetric  plenum  chamber  with  a  flat  circular  faceplate  designed  to  allow  the 
connection  of  either  axisymmetric  or  planar  measurement  sections.  Flow  leaving  the  attached 
measurement  section  is  vented  outside  through  a  large  silencing  duct.  A  Fisher  model  EK-667, 
4  inch  control  valve  equipped  with  a  40%  capacity  cage  is  located  in  the  supply  line.  The 
control  valve,  in  conjunction  with  a  Fisher  type  TL101  electronic  controller,  was  used  to 
regulate  the  plenum  pressure  which  could  be  automatically  held  constant  to  within  ±0.69  kPa 


(±0.1  psi)  as  the  tank  farm  pressure  varied.  A  Wallace  &  Tieman  model  61  A-l A-0200 
precision  absolute  pressure  gage  was  used  to  monitor  the  plenum  chamber  pressure. 

A  small  scale,  planar  two-dimensional  supersonic  wind  tunnel  was  fabricated  out  of 
aluminum  for  use  in  this  investigation.  The  wind  tunnel  consists  of  a  constant  area  stagnation 
chamber  followed  by  a  symmetric  supersonic  nozzle  feeding  directly  into  the  test  section.  A 
symmetric  nozzle  is  used  such  that  the  upper  and  lower  wall  boundary  layers  entering  the  test 
section  have  identical  properties.  A  variable  area  diffuser  is  located  at  the  exit  of  the  test  section 
which  can  be  adjusted  to  yield  either  a  diverging  or  converging-diverging  geometry  .  The  wind 
tunnel  has  a  constant  width  of  7.62  cm  (3.0  in.)  over  the  entire  length.  The  entrance  to  the 
stagnation  chamber  is  mounted  directly  to  the  faceplate  of  the  large  plenum  chamber  and  the 
diffuser  exit  is  connected  to  the  silencing  duct  venting  outside.  Two  interchangeable  nozzle 
blocks  were  designed  and  built,  one  for  Mach  2.45  and  the  other  for  Mach  1.60.  The 
supersonic  nozzles  were  designed  using  the  nozzle  design  program  of  Carroll,  et  al.76  w  hich 
employs  an  inviscid  method  of  characteristics  analysis  to  design  supersonic  nozzles  with 
uniform  exit  flow.  Opdcal  access  to  the  test  section  is  provided  by  removable  windows  in  the 
wind  tunnel  side  walls.  The  40.64  cm  (16  in.)  long  by  5.08  cm  (2  in.)  high,  1.27  cm  (0.5  in.) 
thick  soda  lime  glass  is  mounted  in  a  symmetric  frame  which  can  be  rotated  to  view  either  the 
forw  ard  or  rear  portion  of  the  test  section  with  a  7.62  cm  (3  in.)  overlap  in  the  viewing  areas. 
The  surfaces  of  the  glass  were  not  polished  nor  was  any  flatness  specification  imposed. 
Samples  of  optical  glass  with  a  high  degree  of  flatness  and  parallelitv  displayed  only  slightly 
better  optical  quality  than  the  unpolished  soda  lime  glass  when  viewed  through  the  Schlieren 
system.  For  this  reason,  the  considerably  cheaper  unpolished  soda  lime  glass  was  selected 
As  a  side  note,  heat  treated  Pvrex  glass  was  found  to  have  a  high  level  of  internal  stress  which 
appeared  as  stnations  when  viewed  with  the  Schlieren,  even  when  the  surface  was  highly 
polished  to  a  strict  flatness  tolerance.  A  photograph  of  the  wind  tunnel  is  shown  in  Figure  10 
The  faceplate  of  the  large  axisymmetric  stagnation  chamber  is  visible  to  the  left  and  the 
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silencing  duct  is  on  the  right.  The  near  window  has  been  removed  to  reveal  the  Mach  1 .6  test 
section.  The  far  window  frame  is  mounted  with  the  glass  in  the  forward  position. 

A  description  of  the  wind  tunnel  configuration  with  the  Mach  2.45  nozzle  is  given  first. 
A  photograph  of  the  Mach  2.45  nozzle  blocks  is  shown  in  Figure  1 1.  The  upper  and  lower 
stagnation  chamber,  nozzle,  and  test  section  walls  were  each  fabricated  from  a  single  piece  of 
aluminum  to  avoid  seams  between  sections.  The  stagnation  chamber  had  a  constant  height  of 
10.16  cm  (4.0  in.).  The  nozzle  throat,  located  33.528  cm  (13.2  in.)  from  the  entrance,  had  a 
height  of  1.448  cm  (0.57  in.).  The  nozzle  exit  height  was  3.81  cm  (1.5  in.)  with  the  distance 
from  the  throat  to  the  nozzle  exit  being  8.179  cm  (3.22  in.).  This  nozzle  was  designed  for  a 
Mach  number  of  2.50  but  due  to  boundary  layer  growth  on  the  nozzle  walls  yielded  a  Mach 
number  of  2.45  based  on  wall  static  pressure  measurements.  Following  the  recommendations 
of  Bradshaw77and  Baines  and  Peterson7**,  two  removable  screens  each  with  57.4%  open  area 
(44x44  mesh  stainless  steel  wire  cloth,  0.0055  in.  wire  diameter),  were  located  in  the 
rectangular  stagnation  chamber,  16.36  cm  (6.44  in.)  and  20.19  cm  (7.95  in.)  from  the 
entrance,  to  help  provide  uniform  mean  and  fluctuating  velocity  profiles.  The  screens  are 
mounted  in  aluminum  frames  which  are  slid  into  the  slots  visible  in  the  stagnation  chamber 
portion  of  the  nozzle  blocks  of  Figure  1 1.  Originally,  the  Mach  2.45  test  section  was  made 
with  a  constant  height  of  3.81  cm  (1.5  in.).  However,  initial  Schlieren  visualizations  with  th is 
nozzle  block  revealed  undesirable  oscillatory  movement  of  the  shock  train  location  on  the  order 
of  a  boundary  layer  thickness  in  amplitude.  The  frequency  of  this  shock  motion  was  not 
determined.  In  an  effort  to  eliminate  this  shock  motion,  the  test  section  walls  were  modified 
such  that  the  upper  and  lower  walls  could  be  diverged  by  a  variable  angle  of  up  to  1 .25  degrees 
each  with  the  divergence  beginning  at  the  nozzle  exit.  The  variable  area  exit  diffuser  was  also 
used  to  isolate  the  test  section  from  downstream  fluctuations  by  forming  a  nearly  sonic  second 
throat.  A  series  of  Schlieren  visualizations  with  various  test  section  divergence  angles  and 
diffuser  geometries  showed  that  the  shock  train  unsteadiness  was  not  visibly  affected  by  the 
divergence  angle.  However,  the  second  throat  diffuser  did  tend  to  make  the  mean  shock 
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location  insensitive  to  relatively  large  variations  in  the  plenum  pressure  of  the  order  of  34.5  to 
69.1  kPa  (5  to  10  psia).  For  the  results  presented  in  this  report,  the  test  section  divergence 
angle  was  held  fixed  at  0.25  degrees  for  both  the  top  and  bottom  walls  to  give  a  zero  pressure 
gradient  in  the  fully  started  test  section.  Static  pressure  taps  were  located  along  the  centerline 
of  the  upper  and  lower  test  section  wails.  On  the  upper  wall,  40  taps  were  installed,  the  firsi 
tap  being  at  a  location  0.0762  cm  (0.03  in.)  from  the  nozzle  exit.  A  spacing  between  taps  of 
1.27  cm  (0.5  in.)  was  used  from  0.0762  cm  (0.03  in.)  to  38.18  cm  (15.03  in.)  downstream  of 
the  nozzle  exit,  a  spacing  of  2.54  cm  (1.0  in.)  was  used  from  38.18  cm  ( 15.03  in.)  to  53.42 
cm  (21.03  in.)  downstream,  and  a  spacing  of  5.08  cm  (2.0  in.)  was  used  from  53.42  cm 
(21.03  in.)  to  68.66  cm  (27.03  in.)  downstream.  On  the  lower  test  section  wall,  14  tans, 
spaced  5.08  cm  (2.0  in.)  apart,  were  installed  with  the  first  tap  located  2.62  cm  (1.03  in.) 
downtream  of  the  nozzle  exit.  The  pressure  tap  diameter  on  the  inside  surface  of  the  test  section 
w  as  0.0635  cm  (0.025  in.)  I.  D.  On  the  outside  of  the  test  section  a  0.16  cm  (0.063  in.)  O.  D. 
stainless  steel  tubulation  was  mounted  to  each  tap  opening.  Vinyl  tubing  was  then  used  to 
connect  the  tabulation  to  the  pressure  measuring  device.  Vacuum  grease  was  used  to  seal  all 
mating  components  of  the  wind  tunnel.  Initially,  soapy  water  applied  to  the  joints  of  the  tunnel 
revealed  that  leaks  were  present  around  the  wind  tunnel  window  frames.  Surface  flow  results 
also  showed  the  presence  of  these  leaks.  Additional  mounting  studs  were  added  to  rigidly  hold 
the  window  in  place  and  linear  o-ring  material  was  mounted  in  the  side  of  the  nozzle  blocks. 
These  steps  corrected  the  leakage  problem. 

The  Mach  1.6  nozzle  was  designed  for  a  Mach  number  of  1.615  to  allow  for  boundary 
layer  growth  on  the  nozzle  walls.  Based  on  pressure  measurements,  the  nozzle  exit  Mach 
number  is  \  ?ry  close  to  1.60.  This  nozzle  has  a  throat  height  of  2.54  cm  (1.0  in.)  and  an  exit 
height  of  3.206  cm  (1.262  in.)  with  a  length  of  4.628  cm  (1.822  in.)  from  the  throat  to  the 
nozzle  exit.  The  throat  is  located  37.34  cm  (14.7  in)  from  the  entrance.  The  stagnation 
chamber  has  he  same  height  as  for  the  Mach  2.45  nozzle  and  the  two  flow  screens  are  again 
located  at  the  same  distance  from  the  entrance.  With  the  Mach  1.6  nozzle,  additional  flow 


straightening  in  the  stagnation  chamber  was  provided  through  honeycomb  installed  3.81  cm 
<  1.5  in)  from  the  entrance.  The  honeycomb  had  a  cell  size  of  0.635  cm  (0.25  in.)  and  a  cell 
length  of  6.35  cm  (2.5  in.)  with  the  web  material  being  0.127  mm  (0.005  in.)  thick  stainless 
steel.  The  test  section  extends  from  the  nozzle  exit  to  75.38  cm  (29.67  in)  downstream.  The 
boundary  layer  growth  in  the  Mach  1.6  test  section  was  predicted  using  the  boundary  layer 
code  of  Dutton  and  Addy.79  Based  on  these  results,  a  fixed  divergence  angle  of  0.13  degrees 
was  used  for  both  the  upper  ar.d  lower  test  section  walls  in  an  effort  to  provide  a  zero  pressure 
gradient  in  the  fully  started  test  section.  Similar  boundary  layer  calculations  with  the  Mach 
2.45  test  section  showed  that  this  method  accurately  predicted  the  boundary  layer  growth  for 
that  case.  As  will  be  shown  in  the  results  section,  a  small  pressure  gradient  resulted  in  the 
Mach  1.6  test  section  in  spite  of  these  efforts.  Pressure  taps  were  located  only  along  the 
centerline  of  the  upper  test  section  wall  for  the  Mach  1.6  nozzle.  The  first  tap  was  located  at 
the  nozzle  exit.  The  second  tap  was  located  2.350  cm  (0.925  in.)  downstream  of  the  nozzle 
exit.  From  2.350  cm  (0.925  in.)  to  9.970  cm  (3.925  in.)  the  tap  spacing  was  2.54  cm  (1.0 
in  ).  From  9.970  cm  (3.925  in.)  to  56.960  cm  (22.425  in.)  downstream  the  tap  spacing  was 
1.27  cm  (0.5  in.).  From  56.960  cm  (22.425  in.)  to  72.20  cm  (28.425  in.)  downstream,  the 
tap  spacing  was  2.54  cm  (1.0  in.).  As  for  the  Mach  2.45  nozzle,  the  inside  test  section 
pressure  tap  diameter  is  0.0635  cm  (0.025  in.)  I.  D.  Connection  to  each  tap  is  provided  on  the 
outside  of  the  test  section  by  0.16  cm  (0.063  in.)  O.  D.  stainless  steel  tubulations.  Vacuum 
grease  was  used  to  seal  all  mating  parts  of  the  wind  tunnel.  Soapy  water  applied  to  the  joints 
w  hile  the  tunnel  was  running  did  not  reveal  any  leaks,  with  surface  flow  results  confirming  this 
finding. 

With  both  nozzle  blocks  the  particles  for  the  LDV  measurements  may  be  injected  into 
the  large  axisymmetric  plenum  chamber  or  along  the  vertical  centerline  of  the  rectangular 
stagnation  chamber.  The  vertical  particle  seeder  consists  of  a  0.476  cm  (3/16  in.)  O.D.  tube 
(<*7  hypodermic  tubing)  which  spans  the  10.16  cm  height  of  the  stagnation  chamber  and  is 
located  before  the  two  screens,  and  after  the  honeycomb  for  the  Mach  1.6  nozzle.  Nine 
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aligned,  evenly  spaced  0.132  cm  (0.052  in.)  diameter  holes  are  drilled  in  the  wall  of  the  seeder 
tube.  The  seeder  tube  may  be  rotated  to  point  the  nine  holes  at  any  desired  angle  to  the  mean 
flow.  Seeding  only  the  vertical  centerline  of  the  stagnation  chamber  should  help  prevent 
accumulation  of  seed  particles  on  the  windows.  Complete  information  on  the  type  of  seed 
particles  used  and  the  particle  generator  is  given  in  the  section  on  LDV  procedures. 

B.  SCHLIEREN,  SURFACE  FLOW,  AND  PRESSURE  MEASURING  PROCEDURES 

Schlieren  photography  was  used  to  investigate  qualitative  features  of  the  shock  train 
phenomenon.  A  standard  Toepler  Schlieren  arrangement  was  used.  The  light  source  was  a 
Xenon  model  457  1.4  p.s  duration  light  source  which  could  be  operated  in  either  a  single  spark 
mode  or  in  a  stroboscopic  mode  for  continuous  viewing.  Either  12  inch  or  8  inch  diameter 
parabolic  mirrors  were  used  to  collimate  the  light.  A  box  camera  using  Polaroid  Type  55 
Positive/Negative  4x5  sheet  film  (ASA  50)  was  used  with  the  12  inch  mirror  set.  Both  the  box 
camera  and  a  35  mm  camera  mount  were  available  for  the  8  inch  mirror  set.  A  Nikon  F2 
camera  body  with  Kodak  Pan-X  film  (ASA  32)  was  used  with  the  35  mm  mount.  The  12  inch 
mirrors  were  used  for  the  Mach  2.45  case  to  give  the  maximum  viewing  area.  For  the  Mach 
1.6  case,  the  interaction  was  shorter  and  the  8  inch  mirrors  with  the  more  convenient  35  nr  i 
mount  w’ere  used. 

Surface  flow  visualization  was  used  to  give  information  concerning  the  separation  and 
reattachment  locations  and  the  two-dimensionality  of  the  flow.  A  mixture  of  600  weight  gear 
oil  and  lampblack  was  used.  This  mixture  was  either  spread  evenly  over  the  test  section 
surface  or  placed  in  discrete  dots  on  the  surface.  The  tunnel  was  then  run  for  a  sufficient 
amount  of  time  to  allow  the  surface  streak  pattern  to  set  up.  After  shutting  down  the  tunnel,  a 
sheet  of  tissue  paper  was  carefuly  laid  on  the  surface  to  absorb  the  resulting  pattern.  A 
photocopy  of  the  tissue  paper  was  then  made  to  obtain  a  permanent  record.  While  the  tunnel 
was  running,  notes  concerning  the  direction  of  the  oil  flow  on  the  surface  were  taken. 
Appropriate  arrows  were  latter  superimposed  on  the  pattern  to  show  these  directions.  The  oil 


pattern  was  also  observed  during  the  tunnel  shutdown  process  to  ensure  that  this  transient  did 
not  significantly  disturb  the  pattern. 

Measurements  of  the  wall  static  pressure  were  made  with  two  Pressure  Systems 
Incorporated  (PSD  model  DP  (MOOT  digital  pressure  transmitters.  One  of  these  transmitters  is 
visible  below  the  tunnel  in  Figure  10.  With  this  system  a  separate  transducer  is  connected  to 
each  pressure  port.  The  ports  are  scanned  electronically  and  an  onboard  microprocessor 
reduces  the  data,  giving  gage  pressure,  and  transfers  it  to  a  host  computer.  The  first  DP  (MOOT 
has  sixteen  ±15  psig  and  sixteen  ±30  psig  transducers,  the  second  DP  (MOOT  has  sixteen  ±45 
psig  and  sixteen  ±100  psig  transducers.  These  transducers  were  calibrated  with  a  primary 
standard  and  the  slope  of  each  calibration  curve  was  stored  in  the  transmitter  memory.  A  zero 
point  calibration  was  performed  before  each  scan  to  account  for  varying  atmospheric  pressure. 
The  estimated  accuracy  of  the  two  lower  range  transducers  was  ±  0.0069  kPa  (±  0.01  psi)  and 
the  accuracy  of  the  two  higher  range  transducers  was  ±  0.0172  kPa  (±0.025  psi).  The  two 
pressure  transmitters  were  connected  to  a  HP  9000  host  computer  through  a  manual  switching 
box  such  that  a  scan  with  the  first  transmitter  could  be  immediately  followed  by  a  scan  with  the 
second  transmitter.  A  scan  of  32  ports  took  approximately  5  seconds.  A  complete  scan  using 
both  transmitters  took  a  total  time  of  between  10  to  15  seconds.  All  data  analysis  and  plotting 
was  done  on  the  HP  9000. 

C.  LDV  PROCEDURES 

The  laser  Doppler  Velocimeter  (LDV)  technique  was  selected  as  the  method  to  be  used 
in  making  the  mean  and  fluctuating  velocity  measurements.  Johnson  and  Rose80  and 
Ardonceau81  have  shown  that  both  LDV  and  hot  wire  measurements  have  limitations  in 
transonic  Hows.  Hot  w  ire  measurements  exhibit  large  variations  in  the  mass  flux  sensitivity  in 
the  range  of  Mach  numbers  0.8  <  M  <  1.4.  This  effect  is  caused  by  detached  shock  waves 
from  the  hot  wire  supports  interacting  with  the  hot  wire.  Thus,  large  uncertainties  in  hot  wire 
measurements  occur  at  the  Mach  numbers  present  in  the  shock  train.  Additionally,  the  physical 
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presence  of  the  hot  wire  probe  in  sensitive  transonic  and  separated  flows  may  noticeably 
disrupt  the  flow.  The  LDV  is  seen  as  an  effective  tool  in  the  study  of  separated  shock 
wave/boundary  layer  interactions  provided  the  effects  of  panicle  lag,  velocity  bias,  and  fringe 
bias  can  be  accounted  for.  The  nonintrusive  nature  of  the  LDV  and  the  unambiguous  detection 
of  reversed  flow  ar;  two  strong  motivations  for  using  tne  LDV  technique  in  sensitive  shock 
separated  flows. 

Since  LDV  is  a  relatively  new  and  complicated  technique,  it  will  be  described  in  some 
detail.  The  LDV  system  to  be  used  in  this  study  is  essentially  the  same  system  used  by 
Petne83,  Samimy83,  and  Kuntz84.  The  system  consists  of  a  Spectra  Physics  model  165  argon- 
ion  5  watt  laser,  TSI  optics  and  electroniccs,  a  three-dimensional  traversing  mechanism,  and  a 
PDP  11/73  minicomputer.  The  experience  of  these  previous  researchers  has  contributed  to  a 
valuable  body  of  knowledge  in  applying  the  LDV  to  high  speed  flows.  A  description  of  the 
current  LDV  operating  procedures  follows. 

Dual-beam,  two-component,  coincident  velocity  measurements  are  made  using  two 
laser  beams  of  488  nm  (blue)  and  514.4  nm  (green).  Appropriate  optics  are  used  to  generate  a 
measurements  volume  consisting  of  two  sets  of  orthogonal  fringes  rotated  at  approximately 
±45  deg.  to  the  mean  flow  direction.  The  available  optical  arrangements  are  listed  in  Tables  5 
4,  and  5.  Beam  expansion  is  available  for  beam  spacings  of  22  and  13  mm  prior  to  the 
expander.  The  configuration  selected  for  this  study  is  a  13  mm  beam  spacing  with  beam 
expansion  and  a  350  mm  transmitting  lens.  This  configuration  gives  the  advantages  of  a  small 
measurement  volume  diameter  and  a  relatively  large  beam  spacing.  Frequency  shifting  at  40 
MHz  will  be  used  to  detect  reverse  flow.  Particles  travelling  at  shallow  angles  to  the  fringes 
may  not  cross  the  required  number  of  fringes  for  a  valid  measurement  causing  a  bias  toward 
measuring  the  velocities  of  particles  travelling  more  nearly  perpendicular  to  the  fringes.  This 
bias  is  commonly  called  fringe  bias.  The  fringe  bias  effects  were  investigated  using  the 
analysis  of  Samimy83.  This  type  of  bias  can  be  especially  severe  when  making  two- 
component.  coincident  measurements  since  the  probability  of  a  valid  coincident  measurement  is 
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equal  to  the  product  of  the  probabilites  of  valid  measurements  in  each  channel.  The 
combination  of  frequency  shifting,  large  fringe  spacing  (and  correspondingly  higher  fringe 
velocity  ).  ±45  deg.  fringe  orientation,  and  the  requirement  of  4  fringe  crossings  for  a  valid 
measurement  has  the  destreable  effect  of  virtually  eliminating  fringe  bias.  The  receiving  optics 
are  used  in  the  forward  scatter  mode  at  10  degrees  off  axis.  Off  axis  collection  reduces  the 
effective  length  ot  the  measurement  volume  and  helps  avoid  problems  with  flare  and  stray  laser 
light.  A  pinhole  aperature  and  optical  color  filter  in  front  of  each  photodetector  also  helps 
eliminate  stray  light,  thereby  increasing  the  signal-to-noise  ratio  of  the  system.  Beam 
expansion  also  improves  this  ratio. 

Two  TSI  model  1990B  frequency  counters,  operating  in  the  single  measurement  per 
burst  mode,  are  used  to  process  the  signal  from  the  photodetectors.  The  signal  is  high-  and 
low-pass  filtered  to  eliminate  the  pedestal  and  noise,  respectively,  from  the  signal.  A 
comparison  between  the  time  for  2  and  4  cycles  of  the  PMT  signal  is  used  to  validate  the 
signal,  and  a  master  interface  is  used  to  check  for  coincidence  between  the  two  counter  outputs. 
This  coincidence  window  is  imposed  to  ensure  that  the  measurement  is  made  from  the  same 
panicle  in  each  channel.  Digital  output  from  the  counters  is  transferred  directly  to  the  PDP 
1 1/73  and  saved  on  hard  disk.  Date  reduction  can  be  performed  on  the  PDP  1 1/73  or  the  raw 
data  can  be  transferred  to  a  Hewlett-Packard  9000  series  minicomputer  for  data  reduction  and 
graphical  output.  The  frequency  counters  are  capable  of  sampling  at  a  speed  much  higher  than 
the  data  rate.  Hence,  the  data  is  not  sampled  at  equal  time  intervals.  Rather,  each  data  sample 
which  meets  the  coincidence  requirement  is  collected,  regardless  of  the  time  between  samples. 
As  a  consequence,  the  resulting  data  is  completely  velocity  biased  and  a  two-dimensional 
velocity  magnitude  bias  correction85,  similar  to  that  of  Mclaughlin  and  Tiederman86.  is 
employed  to  correct  the  data.  A  TSI  six  jet  atomizer  is  used  to  generate  seed  particles,  I'sing 
this  atomizer.  Petrie8-  produced  approximately  1  micron  diameter  particles  when  50  cp 


V.  EXPERIMENTAL  RESULTS 


The  results  of  the  integrated  numerical  and  experimental  investigation  of  multiple  shock 
wave/turbulent  boundary  layer  interactions  in  confined  rectangular  ducts  are  presented  in  this 
section.  Two  incoming,  undisturbed  Mach  numbers  were  considered,  Mach  1.6  and  2.45 
For  both  the  Mach  1.6  and  2.45  interactions,  spark  Schlieren  flow  visualization,  surface  oil 
flow  visualization,  and  wall  static  pressure  measurements  w'ere  made.  LDV  measurements  of  a 
Mach  1.6  interaction  are  currently  being  made.  The  Mach  2.45  interaction  was  unsuitable  for 
LDV  measurements  with  the  available  equipment  due  to  shock  train  unsteadiness.  Numerical 
computations  of  the  Mach  1.6  interaction  are  also  being  performed.  The  experimental  results 
are  discussed  in  the  following  sections  and,  when  appropriate,  the  corresponding  numerical 
results  are  also  presented. 

A.  MACH  2.45  RESULTS 

The  experimental  results  for  the  Mach  2.45  shock  train  interaction  are  discussed  in  this 
section.  The  preliminary  experiments  with  the  Mach  2.45  test  section  have  been  previously 
described  in  Chapter  4.  For  the  results  given  here  a  constant  test  section  divergence  angle  of 
0.25°  was  set  on  the  upper  and  lower  walls  giving  a  total  divergence  angle  of  0.50°.  This  was 
done  to  provide  a  neutral  pressure  gradient  in  the  undisturbed  boundary  layer.  Three  cases  are 
considered,  corresponding  to  three  shock  locations  in  the  test  section.  The  experimental 
operating  conditions  are  summarized  in  Table  6.  For  all  three  cases  the  wtind  tunnel  was 
operated  at  a  nominal  stagnation  pressure  of  31 1  kPa  (45  psia)  measured  in  the  large  plenum 
chamber.  The  adjustable  diffuser  was  configured  as  a  converging-diverging  diffuser.  By 
adjusting  the  size  of  the  second  throat,  the  shock  train  was  positioned  at  the  three  desired 
locations.  This  series  of  experiments  is  useful  in  investigating  the  effects  of  flow  confinement, 
as  described  by  the  ratio  6u/h,  with  the  Mach  number  and  unit  Reynolds  number  held  constant. 
The  Reynolds  number  based  on  5U  will  change  slightly  as  5U  changes.  For  these  preliminary 
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results,  the  undisturbed  boundary  layer  thickness,  8U,  defined  as  the  boundary  layer  thickness 
at  the  start  of  the  pressure  rise,  was  estimated  from  the  Schlieren  photograpfs. 


1 .  MACH  2.45  FLOW  VISUALIZATIONS 

The  Schlieren  and  surface  flow'  results  for  the  Mach  2.45  interaction  are  given  in 
Figures  12,  13,  and  14.  In  Figures  12  and  13  the  top  and  bottom  portions  of  the  figures  are 
the  oil  surface  flow-  patterns  of  the  top  and  bottom  test  section  walls,  respectively.  Between 
these  is  shown  the  Schlieren,  to  scale  with  the  surface  flow,  viewed  through  the  window  in  the 
tunnel  side  wall.  As  discussed  below,  the  surface  flow  for  the  third  case  could  not  be 
performed  due  to  flow  unsteadiness.  Therefore,  only  the  side  view  Schlieren  is  given  in  Figure 
14.  In  Figures  13  and  14,  the  Schlieren  picture  shown  consists  of  two  photographs  joined 
together.  The  entire  shock  train  could  not  be  photographed  at  once  due  to  the  limited  field  of 
view  of  the  Schlieren  apparatus.  In  all  three  figures  the  Schlieren  was  set  up  with  a  vertical 
knife  edge  such  that  accelerations  show  as  light  regions  and  decelerations  as  dark  regions.  The 
flow  is  from  left  to  right  in  all  cases. 

A  few  cautionary  words  should  be  made  concerning  the  interpretation  of  these  flow 
visualizations.  The  Schlieren  and  surface  flow  results  are  useful  in  studying  the  qualitative 
features  of  the  flow.  However,  care  must  be  taken  to  not  draw  excessively  quantitative 
conclusions  from  these  results.  In  particular,  caution  must  be  taken  in  judging  the  velocity 
magnitudes  through  the  interaction  based  solely  on  the  Schlieren  photographs.  Also,  the 
surface  flow  results  give  a  good  indication  of  the  flow  directions  at  the  surface,  but  should  not 
be  used  to  evaluate  the  flow  too  far  away  from  the  surface. 

Consider  first  the  case  with  a  confinement  level  of  5u/h  =  0.15,  Figure  12.  In  the 
Schlieren  photograph,  weak  Mach  waves  can  be  seen  in  the  incoming  flow.  These  were 
generated  by  the  growing  boundary  Lver  in  the  nozzle  and  by  any  slight  imperfections  in  the 
nozzle  contour.  The  first  shock  in  the  shock  train  is  an  asymmetric  oblique  shock  pattern.  The 
top  leading  oblique  shock  originates  in  the  upper  wall  boundary  layer  upstream  of  the  location 
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where  the  bottom  leading  oblique  shock  intersects  the  lower  wall  boundary  layer.  The  two 
leading  shocks  cross  below  the  centerline.  After  crossing,  the  two  shocks  terminate  a  fair 
distance  from  the  nozzle  wall  indicating  a  substantial  subsonic  region  near  the  w-all.  Following 
the  first  oblique  shock  pattern,  the  light  region  indicates  that  the  flow  reaccelerates.  It  is  not 
possible  to  determine  if  the  flow  in  this  core  region  is  totally  supersonic  or  if  the  reacceleration 
occurs  from  subsonic  initial  velocities.  The  second  and  later  shocks  are  not  as  well  defined  but 
appear  to  be  nearly  normal  in  character  near  the  lower  wall  and  tend  to  'lie'  along  the  lower 
wall  without  extending  all  the  w-av  to  the  upper  wall.  In  the  shock  pattern  shown,  the  shock 
train  is  described  as  being  'attached'  to  the  lower  wall  even  though  boundary  layer  separation  is 
present,  as  discussed  below.  This  pattern  was  neutrally  stable,  with  the  shock  system 
sometimes  being  attached  to  the  lower  wall  and  sometimes  attached  to  the  upper  wall. 
Occasionally,  the  shock  system  would  flip  from  one  wall  to  the  other  during  the  course  of  a 
run.  A  large  amplification  of  the  size  of  the  turbulence  structure  is  observed  going  through  the 
interaction  and  the  boundary  layer  is  seen  to  thicken  substantially  through  the  interaction. 

Turning  attention  to  the  surface  flow  pattern  for  this  case,  one  observes  two  highly 
different  patterns  on  the  top  and  bottom  walls.  This  is  not  surprising  in  light  of  the  asymmetric 
shock  pattern  described  above.  On  the  top  wall  a  fairly  straight  separation  line  occurs  slightly 
upstream  of  the  separation  on  the  bottom  wall.  Following  the  separation  on  the  top  wall  is  a 
large  three-dimensional  reverse  flow  region.  The  flow  reattaches  near  the  center  of  the  tunnel 
at  two  'star'  patterns  on  either  side  of  the  centerline.  After  these  star  patterns  the  flow  is  seen 
to  head  downstream,  with  some  of  the  flow  closer  to  the  side  walls  turning  and  heading  back 
upstream.  The  flow  reattaches  near  the  side  walls  further  downstream.  This  type  of 
reattachment  could  be  described  as  a  U-shaped  reattachment  with  the  open  part  of  the  U  facing 
downstream.  Substantial  corner  effects  are  seen  on  the  upper  wall  with  the  flow  near  the 
comers  heading  away  from  the  side  walls.  Following  complete  reattachment,  the  flow  is  again 
fairly  two-dimensional.  On  the  bottom  wall  the  separation  line  is  again  nearly  straight.  The 
separation  region  is  much  smaller  and  the  reattachment  line  is  only  slightly  curved.  Again 
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comer  effects  are  seen  but  are  not  large  until  after  reattachment  on  the  bottom  wall  where  the 
comer  vortex  appears  to  grow  at  a  faster  rate  on  the  bottom  wall  with  only  the  middle  third  of 
the  flow  being  two-dimensional  at  the  downstream  edge  of  the  picture.  In  a  study  of 
supersonic  flow  through  a  square  duct  without  shocks,  Davis,  et  al.^7  observed  a  pair  of 
counter  rotating  secondary  flow  cells  centered  about  the  comer  bisector.  Such  a  phenomenon 
could  explain  the  comer  effect  observed  here. 

Moving  to  the  second  case  with  8u/h  =  0.26  (Figure  13),  the  results  are  similar  to  those 
with  5u/h  =  0.15.  However,  the  shock  system  now  tends  to  flip-flop,  from  top  wall 
attachment  to  bottom  wall  attachment,  or  vice  versa,  more  often.  The  case  of  bottom  wall 
shock  train  attachment  is  shown.  In  the  Schlieren  picture  a  small  amount  of  oil  can  be  seen  on 
the  window,  appearing  as  small  streaks  on  the  glass.  The  Schlieren  picture  is  very  similar  to 
the  5u/h  =  0.15  case  of  Figure  12  with  the  asymmetric  oblique  initial  shock,  followed  by  a 
senes  of  shocks  that  are  more  normal  in  character.  The  overall  length  of  the  interaction  has  not 
noticeably  increased.  The  top  wall  surface  flow  again  shows  a  nearly  straight  separation  line 
with  a  U  shaped  reattachment.  However,  the  flow  now  reattaches  closer  to  the  separation  line. 
The  reattachment  near  the  side  wall  is  clearer  in  this  picture,  showing  a  fan  shaped  pattern  at  the 
wall.  The  region  near  the  side  wall  appears  to  be  strongly  influenced  by  the  side  wall  and 
corner  effects.  Thus,  this  fan  shaped  reattachment  near  the  side  wall  may  be  caused  by 
secondary  flow  cells  in  the  comers.  The  flow  following  reattachment  is  again  fairly  two- 
dimensional  on  the  top  wall  with  the  comer  effects  being  small.  On  the  bottom  wall  the 
separation  line  is  nearly  straight  with  reattachment  occurring  earlier  than  on  the  top  wall  and 
being  more  two-dimensional.  The  second  shock  is  seen  to  cause  a  necking  in  of  the  comer 
flow  and  an  accumulation  of  oil  under  the  shock.  The  third  shock  causes  a  similar  but  weaker 
effect.  The  washed  out  region  between  the  first  and  second  shocks  on  the  bottom  wall  is 
presumed  to  be  caused  by  the  reacceleration  of  the  flow  following  the  first  shock.  Efforts  to 
improve  the  recorded  pattern  by  placing  additional  oil  in  this  washed  out  area  and  rerunning  the 
test  w  ere  unsuccessful.  The  oil  quickly  left  this  region  before  the  other  features  of  the  surface 
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flow  could  set  up.  Again  a  large  increase  in  the  corner  effect  is  seen  on  the  bottom  wall 
following  the  interaction. 

At  5u/h  =  0.35  the  shock  train  would  not  stay  attached  to  a  single  wall  long  enough  for 
the  surface  flow  pattern  to  stabilize.  Therefore,  only  the  Schlieren  visualization  was  obtained 
and  is  shown  in  Figure  14.  The  flow  pattern  is  more  symmetric  now,  with  both  the  first  and 
second  shocks  being  oblique  shocks.  The  later  shocks  still  tend  to  attach  to  either  the  upper  or 
lower  w  all  in  a  neutrally  stable  manner.  The  length  of  the  interaction  is  about  the  same  as  in  the 
two  previous  cases. 

In  all  three  cases,  substantial  flow  unsteadiness  was  present.  This  phenomenon  was 
only  studied  visually  with  the  stroboscopic  Schlieren.  Two  main  types  of  unsteadiness  were 
observed.  First  is  the  neutrally  stable  shock  pattern's  tendency  to  flip  from  one  wall  to  the 
other  at  random  times.  This  phenomenon  was  more  pronounced  with  the  shock  near  the  end  of 
the  duct.  This  could  be  due  to  the  increased  level  of  flow  confinement  or  to  a  change  in  the 
characteristics  of  the  subsonic  flow  between  the  shock  train  and  the  converging-diverging 
diffuser.  The  second  type  of  unsteadiness  was  in  the  form  of  higher  frequency  shock  motion 
in  the  flow  direction.  The  magnitude  of  this  motion  was  on  the  order  of  an  undisturbed 
boundary  layer  thickness.  Ikui29  postulated  that  this  type  of  motion  was  caused  by  the 
reasonance  of  the  subsonic  flow  following  the  shock.  The  diffuser  is  believed  to  have  isolated 
the  shock  from  the  downstream  silencing  duct,  such  that  any  reasonance  would  be  due  to  either 
a  transverse  or  longitudinal  type  of  reasonance  in  the  test  section. 

2.  MACH  2.45  PRESSURE  MEASUREMENTS 

Wall  static  pressure  measurements  were  made  for  the  three  Mach  2.45  cases 
summarized  in  Table  6,  with  the  results  presented  in  Figures  15,  16.  17,  and  18.  Figures  15- 
17  present  results  for  pressure  measurements  made  on  the  top  wall  only  with  the  shock 
attached  to  the  bottom  wall.  In  all  these  figures,  the  wall  static  pressure,  p.  has  been 
normalized  by  the  stagnation  pressure  before  the  interaction,  pq.  The  ideal  value  of  p/pn 
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following  a  normal  shock  at  Mach  2.45  is  0.4324.  The  maximum  measured  value  of  p/po  after 
the  shock  train  interaction  is  0.3742,  0.3496,  and  0.3175  for  8u/h  =  0.15,  0.26,  and  0.35, 
respectively.  Thus,  the  pressure  recovery  through  the  shock  train  ranges  from  86.5%  to 
73.4%  of  the  ideal  normal  shock  value.  This  last  value  of  73.4%  may  be  slightly  low  as  the 
pressure  recovery  may  not  have  been  completed  before  the  exit  of  the  duct  in  this  case.  This  is 
seen  in  Figure  15  vs  here  the  abscissa  is  the  actual  downstream  distance  from  the  nozzle  exit;  in 
the  6u/h  =  0.35  case  the  pressure  is  still  rising  at  the  duct  exit.  The  zero  pressure  gradient  in 
the  incoming  boundary  layer  is  evident.  A  comparison  of  the  traces  shows  that  the  initial 
pressure  nse  is  very  similar  for  all  three  cases,  with  the  steep  initial  pressure  rise  followed  by  a 
pressure  plateau.  This  similarity  is  seen  more  clearly  in  Figure  16,  where  the  pressure  traces 
have  been  shifted  such  that  the  downstream  distance  is  measured  from  the  start  of  the 
interaction.  This  agrees  with  the  Findings  of  Waltrup  and  Billig32  in  their  investigation  of 
shock  trains  in  circular  ducts.  This  result  is  also  consistent  with  the  flow  visualizations.  As 
was  noted  in  the  Schlieren  results,  the  structure  of  the  shock  system,  which  causes  the  pressure 
rise,  does  not  change  appreciably  as  5u/h  increases.  The  top  wall  surface  flow  is  seen  to 
change  to  some  extent  but  would  not  be  expected  to  have  a  great  effect  on  the  wall  pressure. 
As  will  be  shown  later,  the  Mach  1.6  wall  pressure  data  collapsed  well  when  the  downstream 
distance  is  measured  from  the  start  of  the  interaction  and  is  normalized  by  the  undisturbed 
boundary  layer  thickness  at  the  start  of  the  interaction,  8U.  The  Mach  2.45  results  have  been 
plotted  in  such  a  manner  in  Figure  17  but  fail  to  collapse  to  a  single  curve.  In  Figure  18  the  top 
and  bottom  wall  pressure  distributions  for  the  5u/h  =  0.26  case  are  shown.  In  this  case  the 
shock  system  is  attached  to  the  bottom  wall.  (The  bottom  wall  measurement  was  actually 
obtained  using  the  pressure  taps  on  the  top  wall  but  with  the  shock  attached  to  the  top  wall. 
This  was  done  to  take  advantage  of  the  Finer  pressure  tap  spacing  on  the  top  wall.)  The  smooth 
pressure  distribution  occurs  on  the  top  wall  where  the  larger  separation  region  is  present.  The 
lagged  pressure  trace  on  the  bottom  wall  was  caused  by  the  normal  shocks  lying  very  close  to 
the  surface.  The  initial  pressure  nse  starts  on  the  top  wall  First  since  the  upper  oblique  shock 
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slightly  leads  the  bottom  oblique  shock.  The  top  and  bottom  wall  pressure  traces  eventually 
coincide  near  the  end  of  the  interaction.  Differences  between  the  two  profiles  occur  only  in  the 
middle  of  the  shock  train  w  here  normal  shocks  lie  along  one  wall. 

B.  MACH  i.6  RESULTS 

Six  cases  are  considered  for  the  Mach  1.6  interaction.  A  summary  of  these  cases  is 
given  in  Table  7.  The  stagnation  pressure  was  held  at  a  nominal  value  of  20.7  kPa  (30  psiai 
for  this  set  of  experiments  and  the  shock  location  was  varied  by  adjusting  the  throat  gap  on  the 
adjustable  diffuser,  configured  as  a  symmetric  converging-diverging  second  throat  diffuser. 
The  test  section  was  built  with  a  constant  divergence  angle  of  0.13°  on  each  side,  yielding  a 
0.26°  total  divergence  angle.  This  divergence  angle  was  expected  to  provide  a  neutral  pressure 
gradient  in  the  undisturbed  boundary  layer  but  a  small  adverse  pressure  gradient  is  present. 
Therefore,  the  Mach  number  is  approximately  constant  in  this  series  of  experiments,  the  unit 
Reynolds  number  is  held  constant,  and  the  flow  confinement  parameter,  5u/h,  is  varied.  The 
undisturbed  boundary  layer  thickness,  8U,  was  estimated  from  Schlieren  photographs  with  a 
horizontal  knife  edge.  An  exact  value  of  5U  will  be  determined  from  the  LDV  measurements 
once  these  results  are  available.  In  the  following  section,  the  Schlieren  results  are  discussed 
first.  Surface  flow  visualization  for  each  case  was  so  similar  that  only  the  surface  flow  results 
for  Su/h  =  0.34  are  presented.  The  wall  static  pressure  measurements  are  then  described. 

1.  MACH  1.6  FLOW  VISUALIZATIONS 

The  Schlieren  results  are  presented  in  Figures  19  and  20  with  the  flow  from  left  to  right 
and  the  knife  edge  vertical  with  decelerations  showing  as  dark  areas  and  accelerations  as  light 
regions.  In  Figure  19  the  Schlieren  was  adjusted  to  a  lower  sensitivity  to  view  the  shape  of  the 
shock  waves  in  the  center  of  the  duct  without  interference  from  the  side  wall  shock 
wave/boundary  layer  interaction.  The  top  picture  shows  the  shock  positioned  near  the  nozzle 
exit  where  the  boundary  layer  is  thinnest.  Proceeding  from  top  to  bottom,  the  shock  is 
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progressively  further  from  the  nozzle  exit  and  the  boundary  layer  is  correspondingly  thicker, 
yielding  increased  values  of  the  confinement  parameter.  In  Figure  20  the  same  sequence  of 
cases  is  given,  now  with  increased  Schlieren  sensitivity.  The  turbulence  structure  is  more 
visible  in  Figure  20.  but  the  shock  structure  is  somewhat  obscured  by  the  side  wall  interaction 
w  hich  appears  as  a  dark  region  in  front  of  the  first  shock  and  behind  the  central  part  of  the  later 
shocks.  The  shock  train  is  very  symmetric  in  all  cases  except  5u/h  =  0.06  where  a  slight 
asymmetry  occurs.  In  all  cases  an  increase  in  the  boundary  layer  thickness  and  an 
amplification  of  the  turbulence  scales  is  observed  going  through  the  interaction.  The  first 
shock  is  always  bifurcated,  w'hile  the  following  secondary'  shocks  are  not.  The  bifurcated 
shock  consists  of  a  leading  oblique  shock,  a  nearly  normal  outer  shock,  and  a  trailing  oblique 
shock  w  ith  the  three  intersecting  at  the  bifurcation  point.  The  outer  normal  shock  is  concave 
facing  upstream.  A  slip  line  is  also  generated  at  the  bifurcation  point  and  extends  downstream. 
This  slip  line  is  visible  in  Figure  16(a)  for  8y/h  =  0.34  where  the  knife  edge  is  horizontal.  A 
weak  reflected  wave  is  caused  by  the  trailing  shock's  intersection  with  the  boundary  layer.  The 
outer,  nearly  normal  shock,  trailing  shock,  and  reflected  wave  form  a  diamond  shaped  region 
where  the  flow  is  accelerating,  as  indicated  by  the  lighter  shading.  Downstream  of  the 
diamond  region,  the  darker  shading  indicates  the  flow  is  only  weakly  accelerating  or  even 
decelerating.  The  size  and  shape  of  the  first  shock  and  the  diamond  region  remains 
approximately  the  same  as  8y/h  increases,  with  the  height  of  the  bifurcation  point  above  the 
wall  increasing  only  slightly.  The  distance  from  the  end  of  the  diamond  (i.  e.  where  the 
reflected  waves  cross)  to  the  second  shock  increases  as  confinement  levels  increase. 

The  shocks  following  the  first  shock,  termed  secondary  shocks,  are  similar  to  each 
other  but  are  different  in  character  than  the  first  shock.  The  secondary  shocks  are  unbifurcated. 
The  outer  region  is  nearly  normal,  being  concave  facing  downstream  as  opposed  to  the  first 
shock  in  which  the  normal  portion  was  concave  facing  upstream.  An  inflection  point  in  the 
'hock  shape  is  seen  to  occur  where  the  slip  line  from  the  bifurcation  point  of  the  first  shock 
crosses  the  secondary  shocks.  This  is  especially  visible  in  Figure  21(a)  with  horizontal  knife 


59 


edge.  Both  high  and  low  sensitivity  Schlieren  photographs  with  a  vertical  knife  edge  are 
shown  in  Figures  21(b)  and  21(c),  respectively,  such  that  an  easy  comparison  between  the 
three  types  of  Schlieren  photographs  can  be  made.  A  diamond  shaped  reacceleration  region 
follows  each  secondary  shock,  but  is  not  as  well  defined  as  the  first  diamond.  Within  a  single 
shock  train,  the  spacing  between  successive  downstream  shocks  decreases.  As  confinement 
increases,  the  number  of  secondary'  shocks,  the  spacing  between  respective  shocks  in  each 
shock  train,  and  the  overall  length  of  the  interaction  increase. 

The  steadiness  of  the  interaction  was  observed  with  the  Schlieren  light  source  in  the 
continuous  viewing  mode.  The  mean  shock  train  location  was  fixed  with  only  very  small 
amplitude  oscillations  in  the  leading  shock's  location.  Each  successive  secondary  shock 
showed  motion  of  increasing  amplitude.  The  last  few  secondary  shocks  in  each  shock  train 
were  very  weak  and  displayed  movement  on  the  order  of  an  undisturbed  boundary  layer 
thickness,  5U.  This  shock  motion  was  in  the  streamwise  direction.  No  noticeable  motion  was 
observed  in  the  direction  normal  to  the  wall.  These  observations  of  the  shock  unsteadiness 
agree  w  ith  those  of  Ikui,  et  al.-^. 

The  surface  How  visualization  for  5u/h  =  0.34  and  Mach  1.6  is  shown  in  Figure  22. 
As  was  done  with  the  Mach  2.45  case,  the  top  and  bottom  test  section  wall  surface  oil  flow 
streak  patterns  are  shown  above  and  below'  the  side  view  Schlieren  photograph.  This  figure  is 
representative  of  all  the  surface  flow  results  at  Mach  1.6.  On  the  top  wall,  discrete  drops  of  oil 
were  placed  on  the  duct  wall.  On  the  bottom  wall,  the  entire  surface  was  coated  with  oil.  A 
small  separation  region  is  seen  below  the  first  shock.  A  distinct  separation  and  reattachment 
line  were  not  formed  due  to  the  small  size  of  the  separation.  A  slightly  larger  separation  is  seen 
near  the  comer  and  along  the  centerline.  The  secondary  shocks  do  not  appear  to  separate  the 
flow  A  corner  effect  is  set  up  afler  the  side  wall  separation  under  the  first  shock  and  continues 
to  grow  with  downstream  distance.  The  flow  after  the  first  shock  appears  to  flow  out  of  the 
comer,  possibly  due  to  some  type  of  rotating  secondary  flow  in  the  comer.  At  the  dow  nstream 
edge  of  the  picture,  the  middle  50T  of  the  flow  is  still  twodimensional. 


MACH  1.6  PRESSURE  MEASUREMENTS 


The  wall  static  pressure  measurements  for  the  six  cases  in  Table  7  are  presented  in 
Figures  23.  24,  and  25.  The  wall  static  pressure,  P,  has  been  normalized  by  the  upstream 
stagnation  pressure,  po.  In  Figure  23  the  wall  pressure  is  plotted  versus  the  distance  from  the 
nozzle  exit.  The  small  pressure  gradient  in  the  undisturbed  boundary  layer  is  apparent.  As  the 
shock  train  is  moved  away  from  the  nozzle  exit  to  a  location  with  higher  confinement,  the 
iength  of  the  pressure  rise  increases  and  the  overall  pressure  recovery  decreases.  The  ideal 
pressure  recovery  through  a  normal  shock  at  Mach  1.6  is  p/pp  =  0.6634.  The  pressure 
recovers  through  the  shock  train  interaction  was  p/po  =  0.6094,  0.5962,  0.5735,  0.5582, 
0.5358.  and  0.4998  or  91.9%,  89.9%,  86.6%,  84.1%,  80.8%,  and  73.35%  of  the  ideal  value 
for  Su/h  =  0.06.  0.14.  0.20.  0.26,  0.34,  and  0.44,  respectively.  The  pressure  recovery  has 
been  completed  by  the  duct  exit  for  the  first  four  cases,  possibly  the  fifth  case,  but  not  for  the 
sixth,  case  where  the  shock  train  is  closest  to  the  exit.  Two  causes  are  responsible  for  the  drop 
m  pressure  recovery  at  higher  confinement  levels.  First,  the  shock  system  is  longer  with 
higher  losses  for  the  increased  confinement  and  second,  the  interaction  occurs  at  a  slightly 
lower  Mach  number  as  the  boundary  layer  thickens  due  to  the  slight  adverse  oressure  gradient 
in  the  undisturbed  boundary  layer.  In  Figure  24  the  pressure  traces  have  been  shifted  with  the 
downstream  distance  now  being  measured  from  the  start  of  the  interaction.  The  pressure  traces 
fail  to  collapse  to  a  single  curve  as  they  did  in  the  Mach  2.45  case,  the  reason  being  that  the 
spacing  between  shocks  and  the  number  of  shocks  in  the  shock  train  is  now  strongly  linked  to 
the  level  of  confinement.  The  differences  in  the  shock  structure,  i.  e.  an  oblique  leading  shock 
for  M(1  =  2.45  and  a  normal  shock  for  Mu  =  1.6,  could  also  be  a  contributing  factor  for  this 
difference  in  the  pressure  traces.  Recall  that  a  pressure  plateau  was  observed  with  the  oblique 
leading  shock  but  is  no  longer  present  with  a  bifurcated  normal  shock.  Replotting  the  data  w  ith 
the  downstream  distance  again  measured  from  the  start  of  the  interaction  and  now  non- 
dimensionalized  by  the  undisturbed  boundary  layer  thickness,  5U  ,  (Figure  25)  causes  the 


pressure  traces  to  collapse  reasonably  well  to  a  single  curve  with  the  exception  of  the  Su/h  = 


0.06  case.  The  reason  for  this  exception  is  unclear  but  could  be  due  to  the  asymmetric  nature 
of  the  shock  train  for  this  case. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


In  this  final  section  conclusions  are  drawn  from  the  experimental  and  numerical  results 
and  some  recommendations  for  further  research  are  made.  It  is  to  be  emphasized  that  the 
results  reported  herein  are  partial  results,  as  the  LDV  investigation  of  the  Mach  1.6  interaction 
and  numerical  investigation  of  the  multiple  shock  interaction  have  not  yet  been  completed. 
Thus,  only  partial  conclusions  may  be  made  pending  the  completion  of  the  investigation. 

Multiple  shock  wave/turbulent  boundary  layer  interactions  at  incoming  Mach  numbers 
of  2.45  and  1.6  were  studied  experimentally.  The  Mach  2.45  interaction  was  found  to  be 
asymmetric  with  the  leading  shock  being  oblique  in  nature  and  the  following  shocks  being 
more  normal  in  character  but  tending  to  exist  near  either  the  upper  or  lower  duct  wall  but  not 
extending  completely  to  the  opposing  wall.  The  interaction  was  neutrally  stable  with  the 
repeated  shock  pattern  tending  to  alternate  between  lying  along  the  upper  and  lower  duct  walls 
at  random  times.  A  small  scale,  high  frequency  oscillation  in  the  shock  location  was  also 
observed.  The  shape  and  extent  of  the  separation  regions  on  the  upper  and  lower  walls  were 
vastly  different.  The  Mach  2.45  interaction  w'as  deemed  unsuitable  for  the  planned  LDV 
investigation  due  to  the  flow  unsteadiness  and  asymmetry.  Therefore,  an  investigation  of  a 
Mach  1 .6  interaction  was  also  initiated. 

The  Mach  1.6  interaction  consisted  of  a  senes  of  symmetric,  nearly  normal  shocks. 
The  initial,  bifurcated  shock  caused  a  small  separation  region  at  its  foot  while  the  weaker, 
unbifurcated,  secondary  shocks  did  not  separate  the  boundary  layer.  At  Mach  1.6  the 
interaction  was  much  steadier  than  for  the  Mach  2.45  case.  A  detailed  LDV  investigation  of  the 
Mach  1.6  interaction  for  one  shock  train  location  is  currently  being  performed.  These  results 
will  reveal  the  exact  nature  of  the  flow  reacceleration  following  each  shock  and  will  also  give 
information  on  the  turbulence  structure  through  the  interaction. 

Numerical  results  have  been  presented  for  three  test  cases.  The  first  two  test  cases 
were  for  equilibrium  Hat  plate  boundary  layers  and  were  used  to  ensure  that  the  computer  axle 
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was  operating  correctly  and  that  the  computational  method  was  capable  of  predicting  these 
relatively  simple  flows.  The  third  test  case  was  a  separated  oblique  shock/turbulent  boundary 
layer  interaction  at  a  compression  comer.  This  case  is  being  used  to  further  validate  the 
computational  code  and  to  gain  experience  in  calculating  a  more  complicated  flow  field  before 
computing  the  shock  train  case.  The  calculation  for  the  third  test  case,  using  the  Baldw In¬ 
formix  turbulence  model,  did  not  accurately  predict  the  size  of  the  separation  region  caused  by 
the  shock.  This  result  was  not  unexpected  and  has  been  previously  observed  and  reported  in 
the  literature.  A  calculation  for  this  case  with  the  Wilcox-Rubesin  two-equation  turbulence 
model  is  currently  being  performed  and  is  expected  to  yield  much  improved  results.  The  final 
calculation  to  be  performed  is  for  the  Mach  1.6  multiple  normal  shock/turbulent  boundary'  layer 
interaction  being  invesdgated  exoenmentally  with  the  LDV  system. 

Several  areas  requiring  further  research  have  been  identified.  The  steadiness  of  both 
the  multiple  normal  and  oblique  shock/turbulent  boundary  layer  interactions  in  nearly  constant 
area  ducts  is  questionable.  The  exact  mechanisms  which  trigger  and  amplify  the  shock  motion 
are  poorly  understood.  The  relative  steadiness  of  planar  versus  axisymmetric  geometries  is 
aoo  unknown.  Three-dimensional  effects  present  in  planar  geometries  need  to  be  investigated 
further.  Finally,  more  work  needs  to  be  done  at  higher  Mach  numbers,  i.e.  Mach  numbers 
above  1.8.  The  bulk  of  the  previous  research  has  focused  on  low  Mach  number  interactions, 
perhaps  because  the  multiple  shock  interaction  is  inherently  steadier  at  lower  Mach  numbers 
and.  consequently,  e  sier  to  handle  experimentally.  However,  many  applications  will  require 
an  understanding  of  the  shock  train  phenomenon  at  higher  Mach  numbers.  The  transition  ol 
the  'hock  train  from  a  normal  shock  train  to  an  oblique  shock  train  at  higher  Mach  numbers  is 
poorly  understood.  Additionally,  the  exact  nature  of  the  oblique  shock  system  (strong  or  weak 
oblique  shocks.)  at  higher  Mach  numbers  is  unclear. 
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=  514.4  nm  (Green) 


TABLE  6 


MACH  2.45  EXPERIMENTAL  CASES 


5U  (mm) 

5u/h 

Mu 

Po  (psia) 

Re/m  (nr1) 

Regu 

3.0 

0.15 

2.45 

44. 5±.  1 

30.  x  106 

9.0  x  104 

5.4 

0.26 

2.45 

44. 9±.  1 

30.  x  106 

16.2  x  104 

7.3 

0.35 

2.45 

44. 9±.  1 

30.  x  106 

21.9  x  104 

TABLE  7 

MACH  1.6  EXPERIMENTAL  CASES 


5U  (mm) 

5u/h 

Mu 

Po  (psia) 

Re/m  (nr1) 

Regu 

0.9 

0.06 

1.60 

29.91.1 

30.  x  106 

2.7  x  104 

2.3 

0.14 

1.59 

29.91.1 

30.  x  106 

6.9  x  104 

3.3 

0.20 

1.57 

29.91.1 

30.  x  106 

9.9  x  104 

4.3* 

0.26 

1.56 

29.91.1 

30.  x  106 

12.9  x  104 

5.7 

0.34 

1.55 

29.91.1 

30.  x  106 

17.1  x  104 

7.5 

0.44 

1.54 

29.91.1 

30.  x  IQ6 

22.5  x  104 
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Normal  Shock 


Figure  3.  Confined,  separated  single  normal  shock/turbulent  boundary 
layer  interaction 


Figure  6.  Turbulence  kinetic  energy  profile,  equilibrium  flat  plate  boundary  layt 
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APPENDIX  A 


Collected  here  are  the  definitions  of  the  4*,  and  Q;  which  are  required  in  the 
compatibility  relations  of  the  reference  plane  characteristic  outflow  boundary  condition.  The 
4/1  are  treated  as  source  terms  and  all  derivatives  contained  within  the  4*,  are  evaluated  using 
finite  difference  expressions  which  maintain  second  order  accuracy.  The  Oj  are  simply 
convenient  groupings  of  terms. 

^1=}  {y^^Pr]  +  Puti>  •  x^vPt!  +  pV  }  (A1) 


^2  =  - }  { u^- y 4u  +  x|v)  -  y^/p } 

+  2p)3T +  X3y]^  ■  ys[a  +  2p)3T  +  ^1, 

r  3v  3  ul  r  3v  3  ul  1 

•  Hp3r+p37Jc  xUp^+p3FJiij 


1*3  *•  J-  {  vn(-y^u  +  +  x5pn/p  } 

+  ^W(X  +  +  X3xk  +  +  2[L)Zy  + 

3v  3  ul  r  3v  3  ul  1 

+  '  yk[^  +  ^^yJn} 

1  a2  r  i 

^4  =  -  7  { (-uy^  +  vx^lp^  j  -  y  {  (-uy^  +  vx^  j 

7r{a  +  2p)(3r?  *a*  2w(37)2  +  4 (37J  +  p 'Mj 
+  2X  (3^37) +  2p  (37J37) +  yi  [po^ '  y5  [p“^l1 
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